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Abstract 


This  work  deals  with  the  problem  of  obtaining  an  open-loop  solution 
to  the  minimum  delay  dynamic  routing  problem. 

The  dynamic  routing  problem  is  stated  using  a  dynamic  model  suggested 
in  previous  works.  This  work  uses  some  previously  known  properties  of  the 
optimal  solution  and  formulates  the  routing  problem  as  a  cubic  optimization 
problem.  In  general  such  problems  are  very  hard  to  solve;  however,  the 
specific  problem  at  hand  is  finally  formulated  as  a  nonconvex  quadratic 
program  by  using  its  special  structure. 

Two  different  approaches,  based  on  the  latter  representation  of  the 
problem,  are  proposed: 

(a)  Utilization  of  existing  methods  for  solving  nonconvex  quadratic 

programs . 

(b)  Development  of  a  special  purpose  algorithm. 

The  algorithm  is  developed  for  single  destination  networks  with  unity  weight¬ 
ings  in  the  cost  functional,  and  it  finds  the  optimal  solution  by  solving 
a  series  of  linear  programs. 

The  algorithm  is  based  on  a  series  of  specially  developed  theorems. 
These  theorems  provide  us  with  new  insight  into  the  behaviour  of  the  dynamic 
routing  in  networks . 

The  method  is  implemented  by  a  computer  program  and  several  examples 
are  run  to  test  its  applicability. 
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Glossary  of  Notations 
set  of  network  nodes 

set  of  network  nodes  not  including  the  destination  node 

set  of  network  links 

directed  link  from  node  i  to  node  k 

capacity  of  link  (i,k) 

vector  of  state  variables 

sum  of  the  elements  of  vector  x 

vector  of  control  variables 

vector  of  control  variables  on  [t.  ,,t.) 

1  i-l  l 

set  of  states  travelling  on  interior  arcs  on  [t  ,t  ) 

p— l  p 

set  of  states  travelling  on  boundary  arcs  on  [t  ^.t^) 

destination  node 

set  of  numbers  {1,2, ...,f} 

switching  time 

time  period  from  t^  j  to  t^ 

vector  of  variables  equal  to  u^t^ 

T 

plane  with  axes  a  x  and  t 

variable  referring  to  the  optimal  solution 

collection  of  nodes  k  such  that  (i,k)  e  L 

collection  of  nodes  l  such  that  (t,i)  e  L 
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Chapter  1 

Introduction 

1. 1  Basic  Concepts 

This  work  considers  the  problem  of  dynamic  routing  in  data-cotamunica- 
tion  networks,  in  the  sense  that  the  routing  depends  on  the  instantaneous 
queue  length  at  the  network  nodes.  A  dynamic  model  for  this  problem  has 
been  proposed  in  [G7]  and  closed-loop  feedback  optimal  solutions  for  this 
model  have  been  investigated  in  [G7],  [G4],  [G5],  [G6]  and  [G2 ] .  It  has 
been  shown  that  the  full  closed-loop  solution  requires  the  construction  of 
polyhedral  convex  cones,  named  "feedback  control  regions",  with  the  property 
that  the  optimal  control  set  is  constant  within  each  region. 

In  this  work  we  present  a  different  approach  whereby  we  seek  the 
open- loop  optimal  control  for  the  dynamic  model  of  [G7].  In  this  case  the 
controls  depend  only  on  the  initial  condition  (i.e.,  the  initial  congestion 
of  the  network) . 

In  this  work  we  will  propose  some  alternatives  for  solving  this 
problem,  all  of  them  based  on  a  change  of  variables  (proposed  by  F.  Moss) 
in  Segal l's  dynamic  model.  (For  other  approaches  to  the  open- loop  problem 
see  [G4],  pp.  45-49). 

1.2  The  Dynamic  Model  of  a  Data  Communication  Network  [G7] 

1.2.1  Topological  Representation 

A  data  communication  network  may  be  represented  as  a  set  of  nodes 
interconnected  by  a  number  of  links  where  to  each  node  are  connected  a 


4 


number  of  "users"  (see  Figure  1.1).  The  links  are  communication  channels 
used  for  the  transmission  of  information  from  node  to  node,  which  are  assumed 


In  a  network  consisting  of  n+1  nodes  we  associate  with  each  node 
an  integer  in  the  set: 

N  =  {1,2, .. . ,n,n+l}  . 

We  now  denote  (see  Figure  1.2): 

1.  (i,k)  is  the  link  connecting  node  i  to  node  k. 

2.  L  =  { (i,k)  such  that  i,keN  and  there  exists  a  link  from  i  to  k) 

that  is,  L  is  the  collection  of  all  links  in  the  network. 

3.  =  (capacity  of  link  (i,k)  in  units  of  traffic/unit  time,  (i,k)cL). 

4.  E(i)  =  (collection  of  nodes  k  such  that  (i,k)  c  L,  VieN)  . 

5.  I(i)  =  (collection  of  nodes  £  such  that  (£,i)cL,  VieN)  . 


The  users  at  each  node  in  the  network  may  input  messages  whose 
destination  is  any  of  the  other  nodes  in  the  network.  We  denote  this  flow 
input  to  the  network  by: 
j  A 

(t)  =  {rate  of  traffic  with  destination  j  arriving  at  node  i 
(from  associated  users)  at  time  t}  . 


A  particular  message  may  pass  through  a  number  of  nodes  before 
arriving  at  the  destination  node.  Since  the  capacity  of  the  links  is  finite, 
every  message  that  arrives  at  a  node  must  wait  in  a  queue  for  being  forwarded. 
Once  a  message  reaches  its  destination  node,  it  is  immediately  forwarded 
to  the  appropriate  user  without  further  storage.  The  principal  source  of 
delay  in  the  transmission  of  information  in  the  network  is  the  queueing 
delay  in  the  nodes. 


For  each  node  we  define  the  following  state  variables: 
i  A 

x^(t)  ■  (amount  of  traffic  at  node  i  at  time  t  whose  final 


destination  is  node  j,  where  i,jeN,  i  +  j)  . 
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Finally,  we  will  define  the  control  variables  of  the  state  space  model  as: 

u^  ■  {flow  of  information  on  the  link  (i,k)  at  tine  t  with 
final  destination  j>  . 


1 . 2. 3  Dyn amic_ Eguat ions _ and  Constraints 

Assuming  that  the  inputs  are  deterministic  functions  of  time,  the 
time  rate  of  change  of  the  number  of  messages  contained  in  each  node  is 
given  by 

xj(t)  *aj(t)-  l  u^.(t)  +  l  uj.  (t)  ,  7  i»j  e  N,  j  t  i  ,  (1.1) 

1  1  kcE(i)  1K  *el(i) 

W 

with  the  control  and  state  constraints: 


xj(t)  >  0  ,  vt  , 


(1.2) 


u?v(t)  -  0  *  :t  • 


(1.3) 


I  uik(t)  -  Cik  ’  Vt 
jcN  >  r 


\  (i,k)  e  L  ,  j  =  i 


(1.4) 


where  is  the  capacity  of  the  link  (i,k) 


1.2.4  Performance  Index 


As  indicated  in  Section  1.2.2,  x(t)  is  the  amount  of  message  resid¬ 
ing  in  some  node  at  time  t,  then  the  quantity 


vf 

/  x(t)dt 
t 


(1.5) 


is  the  total  time  spent  in  this  node  by  the  traffic  passing  through  it  during 
[tQ.tf],  when  x(t  )  *  0.  That  is,  expression  (1.5)  is  the  total  delay  in 
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the  node  experienced  by  the  messages  represented  by  x(t).  Hence,  the  total 
delay  experienced  by  all  the  messages  as  they  travel  through  the  network 
during  [tQ,t^]  *s  given  by 


D  “  /  t  l  xj(t)]dt  ,  i,  j  e  N,  j  i  i  ,  (1.6) 

where  t^  is  defined  as  the  time  at  which  all  the  message  storage  state 
variables  x|  go  to  zero.  Priorities  can  be  accommodated  in  the  cost  func¬ 
tional  (1.6)  by  associating  non-equal  weightings  to  the  appropriate 
state  variables,  so  that  we  have 

tf 

J  a  /  [  [  o]xJ(t)]dt  ,  i, j  e  N,  j  i  i  ,  (1.7) 

*o  ^ 

with  t^  defined  as  above. 


1.3  Research  Overview 

1.3.1  Objective  and  Approach  of  the  Research 

The  goal  of  this  research  is  to  obtain  an  algorithm  for  solving  the 
open- loop  optimal  dynamic  routing  problem  in  data  communication  networks. 

For  a  given  initial  condition,  the  open-loop  solution  will  give  the 
optimal  control  that  takes  the  initial  state  to  zero,  so  that  the  total 
delay  is  minimized.  This  is  in  contradiction  to  the  feedback  (closed-loop) 
solution,  where  the  optimal  control  is  to  be  found  for  every  value  of  the 
state. 

The  reason  for  looking  for  an  open- loop  solution  is  twofold.  First, 
it  turns  out  from  (G4 ]  and  [G2J  that  a  closed-loop  solution  requires  a  large 
number  of  calculations  and  of  memory  space  for  networks  of  reasonable  size. 
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Second,  there  are  still  a  number  of  unsolved  problems  related  to  the  repre¬ 
sentation  of  the  closed-loop  solution,  such  as  determination  of  the  differ¬ 
ent  regions  in  the  state  space,  determination  of  the  region  the  initial 
condition  belongs  to,  and  determination  of  the  time  in  which  the  state 
trajectory  reaches  another  region  for  knowing  when  to  change  the  controls. 

The  approach  used  in  this  research  was  to  make  use  of  known  results 
from  previous  works  for  receiving  a  new  formulation  of  the  dynamic  routing 
problem  and  then  by  means  of  a  change  of  variables  transform  it  into  a 
quadratic  program.  This  quadratic  problem  was  solved  by  using  linear  pro¬ 
grams  after  developing  an  appropriate  theoretical  background.  Likewise, 
other  approaches  for  solving  it,  based  on  known  methods  for  solving  quadratic 
programs  are  proposed. 

1.3.2  Synopsis  of  the  Research 

The  purpose  of  this  section  is  to  provide  a  brief  summary  of  the 
remaining  chapters. 

Chapter  2:  In  this  chapter  the  optimal  control  problem  is  defined 
and  a  number  of  results  relating  to  the  feedback  solution,  as  appear  in  [G4] , 
[G5],  [G6]  and  [G7],  are  presented.  Then,  a  summary  of  the  properties  of 
single  destination  networks  with  unity  weightings  in  the  cost  functional  is 
presented';  - - — — - - - - 

Chapter  3:  In  this  chapter  the  optimal  control  problem  defined  in 
Chapter  2  is  presented  as  a  cubic  optimization  problem  and  after  a  change 
of  variables  as  a  quadratic  program.  Then  it  is  demonstrated  that  the  cost 
function  is  not  quasi-convex,  pseudo-convex,  concave  or  quasi -concave. 

These  properties  are  explained  in  Appendix  A. 


Chapter  4:  In  this  chapter  special  properties  of  single-destination 
networks  with  unity  weightings  in  the  cost  functional  are  used  for  develop¬ 
ing  a  theory  that  enables  us  to  construct  an  algorithm  for  solving  the  open 
loop  optimal  routing  problem  for  this  kind  of  networks  by  using  linear 
programs.  The  proposed  algorithm  is  implemented  on  a  computer  in  Appendix  C. 

Chapter  5;  In  this  chapter  a  number  of  alternative  approaches  for 
solving  the  problem  are  proposed,  all  of  them  based  on  existing  techniques 
for  solving  the  nonconvex  quadratic  programming.  These  techniques  are  over¬ 
viewed  in  Appendix  B. 

Chapter  6:  In  this  last  chapter  the  results  of  the  research  are 
analyzed  and  some  ideas  for  further  research  are  proposed. 

1.3.3  Contributions 

The  contributions  of  this  work  are: 

(a)  For  the  first  time  the  open-loop  dynamic  routing  problem  is  investi¬ 
gated. 

(b)  As  previously  said  in  Section  1.2,  the  variables  of  the  original 
formulation  of  the  optimal  dynamic  routing  problem  are  the  state  and 
control  variables.  In  this  work  we  made  use  of  the  fact  that  the 
optimal  solution  is  of  the  bang-bang  type  in  order  to  formulate  the 
open-loop  routing  problem  as  dependent  on  new  variables:  the  con¬ 
trols  and  the  time  these  controls  are  applied.  This  new  formulation 
is  a  cubic  optimization  problem  with  a  nonconvex  constraint  region. 
This  type  of  problems  is  quite  difficult,  but  the  special  form  of 
the  cost  function  and  of  the  constraint  region  permits  us  to 
formulate  the  problem  as  a  quadratic  program. 


(c)  The  properties  of  this  quadratic  program  are  investigated,  and  it  is 
shown  that  the  cost  function  does  not  have  any  of  the  following 
properties:  convexity,  quasi -convexity,  pseudo-convexity,  concavity 
or  quasi -concavity. 

(d)  Since  the  cost  function  does  not  possess  any  special  property  that 
will  allow  to  solve  the  open- loop  dynamic  routing  problem  easily  by 
existing  methods,  two  basic  approaches  for  its  solution  are  proposed. 
The  first  is  based  on  existing  methods  from  nonconvex  quadratic  pro¬ 
gramming,  and  the  second  consists  of  developing  a  special  algorithm. 

(e)  The  existing  methods  for  solving  nonconvex  quadratic  programs  are 
reviewed,  and  the  methods  that  seem  to  be  most  adequate  for  the 
solution  of  our  problem  are  more  deeply  analyzed.  All  results  are 
condensed  in  an  Appendix  that  may  provide  the  reader  interested  in 
nonconvex  quadratic  programming  with  important  first  information, 

of  existing  methods  and  with  an  extensive  bibliography  on  the  subject 

(f)  The  main  part  of  this  work  consists  of  the  development  of  an  algor¬ 
ithm  that  reduces  the  open-loop  dynamic  routing  problem  for  single 
destination  networks  with  unity  weightings  in  the  cost  functional  to 
a  series  of  linear  programs.  This  algorithm  is  developed  based  on  a 
number  of  theorems  that  provide  the  required  theoretical  background. 

The  algorithm  consists  of  first  finding  the  final  time  (i.e.  the 
time  the  network  empties)  for  the  optimal  solution  and  the  non- 
switch-optimal  -solution  (i.e.  the  optimal  solution  among  all  the 
solutions  having  no  switches) .  From  the  theoretical  background 
follows  that  this  can  be  calculated  by  solving  a  small  linear  pro¬ 
gram.  The  fact  that  the  final  time  is  known  enables  us  to  find  the 
one-switch-optimal -solution  by  solving  another  linear  program.  The 
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algorithm  continues  adding  switches  in  each  of  the  intervals  between 
switches.  Each  new  step  is  solved  by  linear  programming.  The  algor¬ 
ithm  will  stop  when  the  addition  of  a  series  of  new  possible  switches 
does  not  decrease  the  value  of  the  optimal  cost. 

(g)  The  only  limitation  on  the  capability  of  the  algorithm  to  solve  large 
networks  are  the  properties  of  the  linear  program  used.  In  this  work 
we  use  a  linear  program  from  an  existing  library  in  order  to  test  the 
algorithm. 


The  Optimal  Control  Problem'  and  its  Solution 


2.1  Introduction 

This  chapter  deals  with  the  fundamental  laws  governing  the  solution 
of  the  optimal  dynamic  routing  problem  in  data  communication  networks. 

Based  on  the  model  described  in  Section  1.2,  we  present  the  optimal 
control  problem,  the  necessary  and  sufficient  condition  for  its  solution 
and  the  properties  of  the  optimal  control.  Then,  in  Section  2.5  we  present 
special  properties  of  single  destination  networks  with  unity  weightings  in 
the  cost  functional. 

In  this  work  we  deal  only  with  networks  without  input  flow,  that  is: 

a?(t)  =  0  ,  Vi,j,t  . 

All  results  in  this  chapter  are  presented  without  proof.  The  proofs 
appear  in  [G2],  [G4],  [G5]  and  [G6] . 

2.2  The  Optimal  Control  Problem 

- To-facilitate  the  -express  Ion -of-  the-^rrohlem  in  mathematical  form. 

we  define  the  column  vectors  x  and  a  (dimension  n) ,  u  (dimension  m, 
with  m  being  the  number  of  controls  to  be  determined),  and  £  (dimen¬ 
sion  r,  with  r  being  the  number  of  links  in  the  network),  which  are 
respectively  concatenations  of  the  state  variables,  cost  functional  weight 
ings,  control  variables  and  link  capacities.  For  a  given  network  topology 
we  define  the  n*»  incidence  matrix  B  as  follows:  associated  with  every 
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state  variable  x|,  there  is  a  row  bT  of  B  such  that 

bTu  -  -  l  u{.+  l  uj  ,  i.jeN,  j^i.  (2.1) 
keE(i)  Zel (i)  11 

B  is  composed  entirely  of  +l*s,  -l's  and  0's,  and  every  column  of  B 
has  at  most  two  non-zero  elements.  If  a  particular  column  has  only  one 
non-zero  entry  then  it  is  -1. 

In  a  similar  way,  we  define  the  rxm  matrix  D:  associated  with 

T 

every  link  (i,k)  is  a  row  d  of  D  such  that 

dTu<C.k  ,  (2.2) 

represents  the  constraint  (1.4).  The  elements  of  D  are  0's  and  l's  only, 
and  each  column  has  precisely  one  +1. 

Based  on  the  above  notations  we  can  state  the  data  communication 
network  dynamic  message  routing  problem  as  follows: 


Problem  2.1 

Find  the  set  of  controls  u  as  a  function  of  time  (and  state,  in 
the  case  of  closed  loop  solution) 


u(t)  ,  "tTIto,ffl  , 

that  will  bring  the  state  from  a  given  initial  condition  x(to) 
x(tf)  ■  0  and  will  minimize  the  cost  functional 

tf  T 

J  -  /  [«  ?(t)]dt  , 


— (2.sy 

x  to 
-o 


(2.4) 


subject  to  the  state  dynamics 


J£L 
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x(t)  «  Bu(t)  , 


and  constraints  on  the  state  and  control  variables 


x(t)  >  0  ,  "t  e  [tQ,tf]  , 


(2.S) 


C2.6) 


Du  <_  C 
u  >  0 


(2.7) 


□  Problem  2.1 


Note  that  (2.3)  -  (2.7)  define  an  optimal  control  problem  with  linear  con¬ 
straints  and  linear  cost  functional. 


2.3  Necessary  and  Sufficiency  Conditions 

Theorem  2.1:  Necessary  Conditions  (see  [G5],  pp.  13-19) 

Let  the  scalar  functional  h  be  defined  as  follows: 

h[u(t),X(t)]  =  XT(t)x(t)  =  XT(t)Bu(t)  .  (2.8) 

A  necessary  condition  for  the  control  law  u*(*)  ell  to  be  optimal 
for  Problem  2.1  is  that  it  minimize  h  pointwise  in  time,  namely 


Xl(t)Bu*(t)  £  XT(t)Bu(t) 


(2.9) 


u(t)  e  U  t  e  [to>tf]  • 

The  costate  X(t)  is  possibly  a  discontinuous  function  which  satis¬ 
fies  the  following  differential  equation 


-dX (t)  ■  adt  +  dn(t)  ,  te  [t^t^]  , 


(2.10) 


where  componentwise  dn(t)  satisfies  the  following  complementary  slackness 
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condition 


xj  (t)dnj (t)  =  0 


dT1i  (*)  i  0 


Vt  e  [tQ,tf] 

i,j  e  N,  j  +  i 


(2.11) 

(2.12) 


The  terminal  boundary  condition  for  the  costate  differential  equa¬ 


tion  is 


X(t^)  =  v  free  , 


and  the  transversality  condition  is 


X*(tf)x(tf)  =  0  . 


(2.13) 


(2.14) 


Finally,  the  function  h  is  everywhere  continuous,  i.e. 


h[u(t'),X(t"))  =  h[u(t+),Xft+)]  ,  vte[to>tf]  .  (2. 15) 


□  Theorem  2 . 1 


The  costate  trajectories  depend  on  the  value  of  the  corresponding 
state  variables.  When  the  state  x^  >  0,  equation  (2.11)  implies  that 
dn?  =  0;  therefore  by  differentiating  equation  (2.1)  with  respect  to  time 
we  obtain 

-x|(t)  =  aJ.  .  (2.16) 

When  the  state  x?  =  0,  its  respective  costate  is  possible  discontinuous, 
depending  on  the  nature  of  nTi!  At  points  for  which  ~n?  is  absolutely 

j  A  dni(t) 

continuous  we  define  vk  =  — —  and  by  differentiating  (2.10)  with 
respect  to  time  we  have,  based  also  on  (2.11)  and  (2.12),  that 

-x|(t)  =  aj  +  pj(t)  ,  pj (t)  <0  .  (2.17) 

On  the  other  hand,  at  times  when  n?  experience  jumps  of  magnitude  An? 
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(see  (2.10)) 

A>j  *  -Anj  ,  (2. 18) 

and  from  equations  C2.ll)  and  (2.12)  it  is  known  that  those  jumps  are  posi¬ 
tive.  Hence 

AXi>0  .  (2.19) 

Theorem  2.2:  Sufficiency  Conditions  (see  [G5],  pp.  19-20) 

The  necessary  conditions  of  Theorem  2.1  are  sufficient. 

□  Theorem  2 . 2 

We  conclude  from  Theorem  2.2  that  all  trajectories  from  x^  to 
x(t^)  =  0  satisfying  the  necessary  conditions  of  Theorem  2.1  are  optimal. 


2 . 4  Characterization  of  the  Optimal  Control 

From  the  inequality  (2.9)  we  see  that  the  optimal  control  u*(*) 
is  received  by  solving  at  every  time  t  g  [t^t^]  the  linear  program: 

u*(x)  =  ARG  MIN  [>T(t) Bu(x) ]  .  (2.20) 

u(x)eU 


The  minimization  can  be  performed  on  one  link  at  a  time.  Consider 
the  link  (i,k)  and  a  possible  set  of  associated  controls 


i-1  i+1 


n+1 


uik’uik’ ' ’ ’ *uik  ,uik  ’ * * * ,uik 


A  given  control  may  appear  in  one  of  the  following  ways: 

1.  ujk  enters  into  exactly  two  state  equations  (see  (2.5)): 


x?(t)  *  -uj.(t)  -  l  uj  (t)  +  l  uj.  Ct)  , 

1  lk  qeE(i)  iq  tcE(i) 

q^k 

xj(t)  =  uj.(t)  +  l  uj.(t)-  l  uJL(t)  . 

*  lk  tel (k)  qeE(k)  Kq 

Mi 


(2.21a) 


(2.21b) 


.***' 
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2.  u?k  enters  into  exactly  one  state  equation: 


ik(t)  *  -uk.  (t)  -  l  u?  (t)  +  l  u] .  (t)  . 

1  lk  q«E(i)  iq  UI(i)  11 


q/k 


(2.22) 


Hence,  all  controls  on  link  (i,k)  contribute  the  following  terns  to  ABu: 
(Aj(t)-xJ(t))uJk(t)+(A^(t)-A^(t))uJk(t)  *  ...  + 

(Ak"1(t)'Xi  1(t))uik1(t)  +  ak+1(t)'Xi+1(t))uik1(t)  +  •**  + 


,,n+l„,.  .n+1  .  n+1 

(Xk  "  Xi  (t^uik  ’ 


(2.23) 


where  Xk ft)  «  0.  Equations  (1.4),  (2.20)  and  (2.23)  determine  the  optimal 
control  law  at  time  t  on  (i,k) : 


1. 


if 


and 


if 


u‘k  *  cik  and  uik  *  0  •  1  *  1  ■ 


(X*(t)  -  A*(t))  <  (Ajj(t)  -  xj(t))  ,  Vj  t  l  , 


(Ak(t)  -  xj(t))  <  0  . 


uikct)  *uik‘(t)  *  •••  *  u"kCt)  =  Sk 


(2.24) 


(2.25) 


u?k(t)  =  0  ’  v  J  ^  [*»4+1»..  -.n] 


(Xk(t)  -  X*(t))  =  (x£+1(t)- A*+1(t))  =  ...  = 


-  (A*(t)  -  x“(t))  <  (xj(t)  -  xj(t))  ,  ij  t  [l.t+l, 


and 


(xk(t)  -  xj(t))  =  (x£+1(t)  -  X*+1(t))  =  ...  =  (X™(t)  -  A“(t))  <  0 


,£+l 


3.  uik(t)  +  uik  Ct)  +  . . .  +  uik Ct:)  i  Cik  ,  (2 . 26) 

u?kCt)  "  0  »  j  t  , 

if 

(\(t)-X*(t))  =  (x£+1(t)-X*+1(t))  =  ...  = 

*  (x”(t)  -x”(t))  =  o  , 

and 

(*j(t)  -  X*(t))  >  0  ,  •  j  i  [*,i+l,...,m]  . 

From  the  optimal  control  law  above  we  conclude  that: 

Conclusion  2.1:  (see  [G4],  pp.  72-73) 

There  always  exists  a  solution  to  Problem  2.1  with  controls  piece 
wise  constant  in  time,  and  trajectory  with  piecewise  constant  slopes. 

□  Conclusion  2.1 

From  Conclusion  2.1  it  is  clear  that  there  exists  an  optimal 
trajectory  which  can  be  characterized  by  a  finite  number  of  parameters. 

We  now  present  a  compact  set  of  notations  for  specifying  these  parameters 

Let  F  be  the  set  [l,2,...,f],  and  tQ,tj,...,t^  be  the 
switching  times  (see  Figure  2.1),  and  then  denote: 

Definition  2.1 

u.  =  {the  control  vector  on  te  [t.  ,,t.)}  .  ieF  . 

-1  1-1  1 

Definition  2.2 

t.  =  (the  time  period  between  t.  .  and  t.,  i.e.,  t.  =  t.  -t.  , 


i 


Figure  2.1:  Illustration  of  u^  and  . 

Definition  2.3 

Bp  =  {xj/x?(t)  «  0  ,  Vte  ttp.j^p)}  »  VpeF- 

Hence,  Bp  is  the  set  of  state  variables  that  take  value  zero  during  the 

time  interval  [t  , ,t  ). 

P~1  P 

Definition  2.4 

B(x)  =  {B1,B2,...,Bf}  . 

That  is,  B(x)  is  the  sequence  of  the  sets  Bp.  B(x)  is  called  the 
boundary  sequence. 


Definition  2.5 


I  =  (xj/xjct)  >  0  ,  Vtc  [t  ,t  )}  ,  peF 


Hence,  Ip  is  the  set  of  state  variables  that  are  not  at  zero  value  during 


the  time  interval  [t  ,  ,t  ). 


2.5  Single  Destination  Networks  with  Unity  Weightings 
in  the  Cost  Functional 

2.5.1  The  Optimal  Control  Problem 

Let  us  denote  the  destination  node  as  the  n+1  node  from  the  set  of 
nodes  N  and  then  define  the  set  of  nodes 

N  =  {1,2,. ..,n}  , 

then,  N  is  the  set  of  all  the  nodes  in  the  network  besides  the  destina¬ 
tion  node.  Since  there  is  only  one  destination  node,  it  is  possible  to 
omit  the  index  that  represents  the  destination  node  in  the  control  and 
state  variables  (see  Sections  1.2  and  2.2).  Therefore  we  will  write  x^ 

instead  of  x?+*,  and  u.,  instead  of  u?.+  *.  Likewise,  in  each  link 
l  lk  lk 

(i,k)  there  is  flow  of  information  towards  a  single  destination  node; 
then,  the  D  matrix  (see  (2.7))  is  the  identity  matrix. 

By  using  the  above  notation  it  is  possible  to  state  Problem  2.1 
for  single  destination  networks  with  unity  weightings  in  the  cost  functional 
(a  =  1)  as  follows: 

Problem  2.2: 

Find  the  set  of  controls  u  as  a  function  of  time  (and  state  in 
the  case  of  closed- loop  solution) 

u(t)  ,  t  e  [tQ,tf]  ,  (2.27) 

that  will  bring  the  state  from  a  given  initial  condition  x(t  )  =  x  to 

-o-o 

x(tf)  *  0  and  will  minimize  the  cost  functional 
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subject  to  the  state  dynamics 

x(t)  =  Bu(t)  , 

and  constraints  on  the  state  and  control  variables 
x(t)  >  0  ,  t  e  [tQ,tf]  , 

U{0  <  u. .  <  C.. }  ,  (i,  j)  e  L 

□  Problem  2.2 


2.S.2  Special  Properties 

Proposition  2.1  (see  [G4],  pp.  244-259,  and  [G2],  pp.  77-83): 

The  costates  X^(t)  are  continuous  functions  of  time  and 
tion  X  (t)  >0  ,  i  e  N,  t. 

□  Proposition  2.1 

Proposition  2.2  (see  [G4],  pp.  263-267,  and  [G2],  pp.  77-83): 

There  always  exists  an  optimal  solution  satisfying  x(t)  • 

□  Proposition  2.2 

Proposition  2.3  (see  [G2],  pp.  51-53): 

The  set  of  all  solutions  to  the  problem: 

min  I  X. x.  ,  X.  >  0  ,  x.  e  I  , 

„  “T  1  1  1  l  p  * 


(2.29) 

(2.30) 
(?.31) 


in  addi- 

.  0  ,  Vt 


satisfying 
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and 


x  *  Bu  ,  u  e  U  , 

*i  i  0  *  xi  6  \  ’ 


x.  =  0  ,  V  x.  e  B 

l  i  P 


is  a  subset  of  the  collection  of  all  optimal  solutions  of  the  same  problem, 

with  {X.  ■  1,  x.  e I  }  . 
i  i  P 


□  Proposition  2.3 


2.5.3  Geometrical  Characterization  of  the  Feedback  Space 

The  feedback  solution  to  the  optimal  control  problem,  as  it  appears 
in  [G4],  [G5],  [G6]  and  [G2] ,  is  based  on  the  construction  of  regions  in 
the  state  space  such  that  to  each  region  corresponds  an  optimal  control 
vector.  The  set  of  all  these  regions  covers  the  state  space.  Therefore 
the  corresponding  optimal  controls  are  the  feedback  solution. 

In  [G4]  an  algorithm  is  presented  for  construction  of  regions  with 
the  property  that  to  all  initial  conditions  lying  in  a  given  region  cor¬ 
respond  the  same  set  of  optimal  controls  and  the  same  boundary  sequence. 

The  geometrical  characterization  of  these  regions  is  a  result  of  the  follow¬ 
ing  theorem: 

Theorem  2.3  (see  [G4],  pp.  114-115): 

The  feedback-control-regions  are  convex  polyhedral  cones  in  11°. 

□  Theorem  2 . 3 

The  algorithm  constructs  the  feedback-control -regions  backwards  in 
time,  where  in  each  stage  it  allows  a  set  of  nodes  to  "leave  the  boundary" 
(that  is,  a  set  of  nodes  is  allowed  to  be  different  from  zero  from  a  given 
time  until  -»,  and  the  other  nodes  remain  at  zero  level  all  the  time). 

(See  [G4] ,  pp.  115-135.) 


In  single  destination  networks  with  unity  weightings  in  the  cost 
functional  there  are  properties  that  facilitate  the  construction  of  the 
regions.  Some  of  these  properties  are: 

Property  2.1  (see  [G4] ,  pp.  126-132,  pp.  2S9-260,  and  [G2],  pp.  84-86) 

For  single  destination  networks  with  unity  weightings  in  the  cost 
functional,  there  are  no  break  walls  in  the  state  space.  That  is,  the 
solution  of  one  stage  of  the  algorithm  for  a  given  set  of  nodes  that  leave 
the  boundary  backwards  in  time  is  constant  in  the  entire  time  interval 
(t  ,-»).  Where  t  is  the  time  when  the  set  of  nodes  leave  the  boundary. 

p  p 

□  Property  2 . 1 

Property  2.2 

We  conclude  immediately  from  Property  2.1  that  if  we  look  at  an 
optimal  trajectory  forward  in  time,  there  are  switches  only  when  a  set  of 
nodes  reach  the  boundary  (i.e.,  a  set  of  state  variables  become  zero). 


□  Property  2 . 2 


Representation  of  the  Optimal  Control  Problem 
as  a  Nonconvex  Quadratic  Program 

3.1  Introduction 

In  the  previous  chapter  it  was  proved  that  there  is  always  an 
optimal  solution  Lo  Problem  2.1  with  piecewise  constant  controls  and  piece- 
wise  constant  slopes  in  the  state  trajectory;  therefore,  this  trajectory 
may  be  represented  by  a  finite  number  of  parameters  (see  Conclusion  2.1). 

In  this  chapter.  Problem  2.1  will  be  formulated  in  such  a  way  that 
the  optimization  will  be  done  over  this  finite  set  of  parameters  with  the 
optimal  solution  fulfilling  Proposition  2.2.  This  formulation  (Problem  3.1) 
is  a  cubic  static  optimization  problem.  Solving  such  problems  turns  out 
to  be  quite  difficult,  but  in  Section  3.3  we  simplify  this  problem  by  reduc¬ 
ing  it  to  a  quadratic  program  (Problem  3.2).  By  means  of  this  representa¬ 
tion  of  the  problem  we  shall  try,  in  the  following  chapters,  to  solve  the 
open  loop  dynamic  routing  problem. 


3.2  Representation  of  the  Optimal  Control  Problem 

by  Means  of  the  Controls  and  the  Time  Between  Switches 

3.2.1  The  Optimal  Control  Problem 

In  Chapter  2  it  was  proved  that  there  is  always  an  optimal  solution 
to  Problem  2.1  with  piecewise  constant  controls  and  with  piecewise  constant 
slopes  in  the  state  trajectory  (see  Conclusion  2.1). 


k 
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Based  on  this  fact  we  may  formulate  Problem  2.1  as  a  static  optimiza¬ 
tion  problem,  with  the  controls  u^,  i  e  F,  and  the  time  periods  x^  ,  ieF 

as  the  optimization  variables.  Solution  of  this  optimization  problem  will 

provide  an  optimal  solution  satisfying  Conclusion  2.1. 

Performance  Index 

From  (2.4)  we  see  that  it  is  possible  to  represent  the  performance 

T  T 

index  J  as  the  area  under  the  curve  a  x(t)  in  the  plane  with  axes  o  x 

and  t  (from  now  on,  this  plane  will  be  denoted  by  2). 

We  will  also  define: 

Definition  3.1 

J^.  =  (The  cost  functional  in  the  time  interval  ^,t^),  ieF}  . 

As  mentioned  above,  the  state  trajectories  have  piecewise  constant 

T 

slopes;  then,  since  a  is  constant,  a  x(t)  also  has  piecewise  constant 
slopes  (see  Figure  3.1). 


Figure  3.1:  Trajectory  with  piecewise  constant  slope,  in  the  Z-plane. 


t 
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T  ^ 


n 


The  slopes  in  each  time  interval  are 


x(t  )  -  x(t  .) 
-  P  -  P-1 

T 

P 


then,  from  (2.5)  and  Proposition  2.1 


p  e  F  , 


(3.1) 


x(tp)  *  x(tp  l)  -  BupTp  ,  peF  .  (3.2) 

We  conclude  from  (3.2)  that  the  knowledge  of  x  ,  t  and  u  for 

-op  -p 

pe  F  is  enough  for  knowing  x(t),  t.  Likewise,  (3.2)  tells  us  that  the 
values  of  the  state  variables  at  the  switching  times  are: 


x  (t^)  =  x(tQ)  +  Bu1T1  , 

x(t2)  =  x (t x )  +  Bu2t2  =  x(tQ)  +  BUjTj  +  Bu2t2  , 


x(tf) 


=  x(to)  +  BUjTj  + . +  BufTf  , 


that  is  to  say 

P 

x(t  )  =  x  +  y  Bu.t.  ,  peF  .  (3.3) 

-  p  -o  .  ,  -1  1 
v  i=l 

The  performance  index  J  will  be  calculated  by  calculating  Jp 
for  each  peF  (see  Figure  3.1).  Note  that  Jp  is  equal  to  the  area  of 
the  trapezoid  in  Figure  3.2,  therefore 


J 

P 


=  f  <?Tx(t)dt  =  +  f(tp_i^Tp 

Yl 


» 


peF.  (3.4) 


•TittP_i) 


«<y 


t  ,  T  t 

P-1  p  p 


Figure  3.2:  as  the  area  of  a  trapezoid* 


From  (3.3)  and  (3.4)  we  conclude  that: 


55  +  l  B^T<  +  *  +  l  B«4TJT«  *  PeF  »  C3.5) 


p  2-  L-o  > 


■1  1  -o 


-1  1  p 


or,  after  algebraic  manipulations 


J  =  ax  t  + 
p  --op 


?TUi^riTp  *  5%tp}  •  p ' F  ■  <3- 


and  the  performance  index  J  is  (see  Figure  3.1) 

r  \  T  .  T^r1.  .  ] 


J  =  7  ^a*x  t  +  a 1  l  Bu.t.t  +  =Bu  t2 

p=i  -  -°  p  -  i-i  -1  1  p  2  'P  P 


(3.7) 


End  Condition 


Based  on  (3.3),  the  end  condition  x(t^)  ®  0  can  be  written  as: 


A. 

*  7  Bu.r.  *  0  . 

-1  i 


(3.8) 


State  and  Control  Constraints 

The  positivity  of  the  state  variables  is  assured  by  the  following 
proposition: 
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Proposition  3.1 

In  the  trajectories  satisfying  Conclusion  2.1,  x(t)  0  for 
te[tQ,tf]  iff  x(tp)  >_  0,  peF. 

Proof 

If  x(tp)  ^0  for  peF  then  x(t)  ^0  for  t  e  [tQ,t^]  (since 
the  slopes  of  the  state  trajectory  are  constant  over  the  interval 
[t  ,,t  ),  peF).  On  the  other  hand,  it  is  obvious  that  if  x(t)  >  0, 

p-1  p  - 

Vt  <  t^,  then  x(t  )  0,  peF. 

□  Proposition  3.1 


From  Proposition  3.1  and  (3.3)  follows  that  it  is  possible  to  write 


the  constraint  x(t)  >_  0  as 

P 

IBVi-'?o’  P  *  l,2,...,f-l 


(3.9) 


i=l 


The  control  constraints  (2.7)  remain  without  change  for  all  the 
controls  u. ,  ieF. 

-l 

In  addition  to  these  constraints,  it  is  required  that 

Ti  1.  0  »  1  e  F  ,  (3.10) 

for  assuring  that  t^  >_  t^  *£-2  —  to  ' 

From  (3.8)  and  (3.9)  it  is  obvious  that  for  each  time  interval 

[ti  ^,t^),  ieF,  there  is  a  change  of  Bu^t^  in  the  state  vector  value; 

therefore  in  [t.  , ,t, )  we  have  x  *  Bu. .  Then  the  differential  con- 
l-l  i  -  -i 

straints  (2.5)  hold. 


Based  on  the  above,  we  can  state  the  optimal  dynamic  routing  prob¬ 
lem  as  the  following  static  optimization  problem: 


Problem  3.1 


Find 


u.  and  t.  ,  ieF  , 
-l  1  ' 


C3.ll) 


that  minimize  the  cost  function 


f 

J-  I 

pal 


T  .  T 
a  x  t_  +  a 
--op  - 


Pf1  1  2 

)  Bu.t.t  +  =Bu  r 

iti  -1  1  P  2  “P  P 


subject  to  the  constraints: 


with  the  amount  of  switches  f  -  1  free. 


(3.12) 


u.  >0,  t  .  >  0 ,  ieF  , 

-i  —  l  — 

(3.13) 

Du.  <  C  ,  ieF  , 

-l  —  - 

(3.14) 

P 

IBuiTi-‘Xn  *  p  =  l,2,...,f-l  , 

(3.15) 

f 

l  Bu.t.  =  -x  , 
i=l  -i  1  -o 

(3.16) 

□  Problem  3.1 


3.2.2  Properties  of  the  Proposed  Problem 

As  mentioned  above.  Problem  3.1  is  a  static  formulation  of  Prob¬ 
lem  2.1.  In  light  of  this  fact,  it  is  possible  to  solve  the  optimal 
dynamic  routing  problem  by  using  known  methods  from  static  optimization. 

There  are  general  and  efficient  methods  for  special  problems  only, 
but  there  are  no  general  methods  for  solving  any  nonlinear  optimization 
problem;  i.e.,  for  each  nonlinear  optimization  problem  we  must  find  the 
adequate  method. 
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The  cost  function  of  Problem  3.1  is  cubic  (see  (3.12))  and  its  con¬ 
straint  region  is  quadratic  (see  (3.15)  and  (3.16)),  i.e..  Problem  3.1  is  a 
cubic  optimization  problem.  From  (3.16)  we  see  immediately  that  the  con¬ 
straint  region  is  nonconvex  since  Definition  A.l  does  not  hold  (see 
Example  3.1). 

A  problem  with  these  properties  is  very  hard  to  solve.  In  the  next 

section  we  will  transform  Problem  3.1  to  a  quadratic  one  by  changing  the 

variables  u.  and  r.  to  new  variables. 

-l  i 

Example  3.1 

Refer  to  the  network  of  Figure  3.3.  It  is  obvious  that  there  are  no 
switches  in  the  optimal  solution,  then  constraint  (3.16)  can  be  written  as: 

Ut  *  -X 

O 

The  constraint  region  is  convex  if  for  u, t,  ■  -x  and  u.t_  =  -x 

1  1  o  i  i  o 

then  [Auj  +  (1-X)u2]  [At1  +  (1-A)t2]  =  -xq  ,  VXe[0,l].  But 

2  2  2 

[AUj  '♦  (1-X)u2)  [ATj  +  (1-X)t2]  =  -Xq  +  [-2Au2t2  +  Xu^  -  X  *  Au^  -  X  u^  -  2X--x-<51 
and  obviously  this  expression  is  not  equal  to  -xq  for  all  Xe  [0.1]. 

x(t)  =  x 

O  0 

O— — O 

Figure  3.3:  Network  for  Example  3.1 

□  Example  3.1 

3. 3  The  Optimal  Control  Problem  as  a  Quadratic  Program 

3.3.1  The  Optimal  Control  Problem 

It  can  be  easily  noted  that  in  all  the  quadratic  and  cubic  elements 
of  Problem  3.1  there  is  a  term  of  the  type  u^.  Therefore  it  is  possible 
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to  transform  Problem  3.1  into  a  problem  with  a  lower  degree  by  defining  the 
following  new  variables: 


Definition  3.2 


z.  =  u.x.  ,  i e  F  . 

-l  -li 

(In  the  next  section  there  is  an  explanation  about  the  meaning  of  z^.) 

Based  on  the  definition  above  the  cost  function  (3.12)  can  be 
written  as: 


J  =  L 


f 

r 

L 

p=l 


T  T 

a  x  t  +  a 
--op  - 


V  1 

)Bz.  T  +  =Bz  T 
-1  p  2  -p  p 


(3.17) 


Likewise,  the  constraints  (3.15)  and  (3.16)  can  be  written  as 
P 

l  Bz.  >  “x0  »  P  =  l,2,...,f-l  , 

i=l  1  ° 

and 


(3.18) 


y  Bz.  =  -x  , 
i=l  -1  -° 


(3.19) 


respectively. 

The  constraints  (3.13)  become 

zi  ^  0,  >_  0,  i  e  F  ,  (3.20) 

since,  if  (3.20)  holds,  from  Definition  3.2  u^  >_  0.  On  the  other  hand  if 
(3.13)  holds,  then  z^  >_  0. 

If  we  replace  u^  from  Definition  3.2  in  (3.14),  the  capacity 
constraints  become 

DZj  -  Ct\  <_  0  ,  i  e  F  .  (3.21) 

In  light  of  these  facts,  the  optimal  dynamic  routing  problem  may  be 


formulated  as: 


-  32  - 


Problem  3.2 


Find 


z .  and  t .  ,  i  e  F 

-l  l 


(with  the  amount  of  switches  f  -  1  not  known)  that  minimize 


f 

d  =  I 

p»l 


P-1 

B...  .  •  . 

=1  -i  p  2  -p  p 


T  T|rr*  1 

a  x  t  +  a  >  Bz.  t  +  -Bz  T 

-  “O  P  -  I  “  -A  t»  O  • 


subject  to  the  constraints 

Zi  10,  T.  >  0,  i  e  F  , 


(3.22) 


(3.23) 


(3.24) 


Dz.  -  CT.  <  0  ,  i c  F  , 

-l  -  l  —  *  * 

V 

1  Bz.  >_  -x  ,  p  =  1,2,..., f-1  , 

i=l  “° 


(3.25) 

(3.26) 


□  Problem  3.2 


(3.27) 


3.2.2  Mea.-;ng  of  the  z  -variables 

- 2 - -p - 

If  we  replace  UpTp  Zp*  expression  (3.2)  becomes 

?<y  ■  ?<y  i’  * 


(3.28) 


i.e.,  BZp  is  the  value  of  the  change  in  the  state  variable  vector  x 
during  t  .  If  we  look  only  at  one  element  of  the  vector  z  (denoted  by 

<*ilA  P 


-  cu{k)p  ‘  xp  • 


(3.29) 


r)  4, 


Then  (z'^r  is  the  value  of  the  change  in  and  (on  opposite 

directions)  because  of  the  information  flow  from  node  i  to  node  k  destin- 

,3  ^ 


ated  to  node  j,  during  (denoted  by  (u^)  ). 
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3.3.3  Progerties_of_Problem_3i2 

The  optimal  dynamic  routing  problem  is  formulated  as  a  static  opti¬ 
mization  problem  in  Problem  3.2.  Its  cost  function  is  quadratic  (see 
(3.23)),  and  the  constraint  region  is  linear  (see  (3.24)-  (3.27));  i.e.. 

Problem  3.2  is  a  quadratic  program. 

Much  is  known  about  quadratic  programming  in  comparison  with  other 
non-linear  optimization  problems,  mainly  when  the  cost  function  is  convex, 
concave,  quasi-convex,  pseudo-convex  or  quasi-concave  (the  definition  of 
these  terms  and  an  explanation  of  their  importance  appear  in  Appendix  A). 

As  a  result  of  that,  we  must  determine  whether  or  not  one  of  these 
properties  holds. 

From  Proposition  A.l  it  is  clear  that  if  a  specific  function  is  not 
quasi-convex,  then  it  is  not  pseudo-convex  or  convex.  Hence,  we  will  show 
that  these  three  properties  do  not  hold  by  showing  an  example  with  a  non- 
quasi-convex  function.  This  example  is  based  on  the  following  network: 

Network  3.1  (see  [G4],  pp.  93-104): 


In  Network  3.1,  as  it  appears  in  Figure  3.4,  there  is  only  one 
destination  node  (node  3) .  Then  it  is  possible  to  omit  the  index  that 
represents  the  destination  node  in  the  controls  and  the  state  variables. 
In  light  of  these  facts,  the  upper  index  will  denote  the  time  interval  in 
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question;  for  example  u*2  is  the  control  uJ2  during  the  time  interval 
[tQ,t1)  (we  will  use  this  notation  henceforth  when  referring  to  single 
destination  networks). 


From  Figure  3.4  it  is  clear  that  the  differential  equations  (2.5) 
for  Network  3.1  are 


Xj(t)  *  -ul3(t) -u12(t) +u21(t) 

f 

x2(t)  =  -u23(t)  +  ul2(t)  "  “gift) 

likewise,  if  the  cost  function  weighting  constants  are  all 
turns  to 


(3.30) 


1,  (3.23) 


•  H 

P=it 


Cxol+Xo2)Tp  + 


P-1 
i=l 


23"' Tp  +  2("Z13‘Z23)Tp 


}• 


(3.31) 


with  x  .  being  the  i-th  element  of  the  initial  condition  vector,  and 
oi  6 

ZJk  being  zik  during  the  time  interval  [t^^t^) . 

□  Network  3.1 


Example  3.2 

T 

Referring  to  Network  3.1  with  the  initial  condition  xq  =  (8,2)  , 
we  will  denote 

z  =  (zP  ZP  zP  tP  1 
-p  lz12*  21,Z13,I23J  ‘ 

It  is  easy  to  verify  that  the  solution  (denoted  by  w*) 

* 

Zj  -  (2,0, 3, 3)  ,  Tj  =  9 

-1 

z2  =  (1,0, 2, 2)  ,  t2  =  16 
is  feasible,  and  that  J(w^)  =  95. 


On  the  other  hand,  the  solution  w 


Zj  *  (0,0, 2,1)  ,  Tj  -  8 


z2  =  (2,0,4, 3)  ,  x2  =  8 


is  also  feasible  and  for  it  the  cost  function  is  J(w  )  =  96. 
Referring  to  those  two  solutions 

J(Xw2 +  (l-X)w2)  =  98.875  for  X  =  1/2  , 


that  is  to  say 


J  (Xw1 +  (l-X)w2)  £  max(J (w1) ,J (w2) }  . 


We  then  conclude  that  Definition  A. 5  does  not  hold  for  all  Xe  [0,1]. 

□  Example  3.2 

In  a  similar  way  it  is  possible  to  show  that  in  general  the  cost 
function  is  not  quasi-concave  (see  Example  3.3). 


Referring  tr  Network  3.1  with  the  initial  condition  xq  »  (4,1)  , 
it  is  easy  to  verify  that  the  solution  (denoted  by  w^) 


Zj  -  (0,0,0. 5, 0.5)  ,  ij  =  5 
l2  =  (1.5, 0,2, 2)  ,  r2  =  5 


is  feasible  and  for  it  the  cost  function  is  J(w  )  «  32.5, 


On  the  other  hand,  the  solution  w 


2,  =  (0.4, 0,0. 6, 0.5)  ,  t,  *  5.595 
-1  1  2 

■  w 

z  *  (0,0, 3, 0.9)  ,  t  *>  3.905 

-z  z  , 

2 

is  also  feasible  and  for  it  the  cost  function  is  J(w  )  ■  32.5125. 
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Based  on  these  two  solutions 

J(Xw 1  *  (l-X)w2)  =  32.49  for  X  =  3/4  , 

that  is  to  say 

JCXw1  +  (l-X)w2)  £  minUCw^.JCw2)}  . 

We  then  conclude  that  Definition  A. 8  does  not  hold  for  all  Xe  [0,1]. 

□  Example  3.3 

It  is  possible  to  resume  this  section  with  the  following  proposition 
Proposition  3.2 

Problem  3.2  is  a  nonconvex  quadratic  program,  with  a  cost  function 
that  is  not  concave,  quasi-convex,  quasi-concave  or  pseudo-convex  even  for 
the  case  of  single  destination  networks  with  unity  weightings  in  the  cost 
function. 

□  Proposition  3.2 

In  Appendix  B  there  is  a  general  review  of  the  methods  for  solving 
nonconvex  quadratic  programs,  and  in  Chapter  5  we  give  propositions  for 
solving  Problem  3.2  by  using  these  methods. 

In  the  next  chapter  we  will  solve  Problem  3.2  for  single  destina¬ 
tion  networks  with  unity  weightings  in  the  cost  function  by  using  special 
properties  of  these  networks. 


3 


Chapter  4 

Open- Loop  Solutions 

for  the  Dynamic  Routing  Problem  for  Single  Destination  Networks 
With  Unity  Weightings  in  the  Cost  Functional 
by  Using  Linear  Programming 


4 , 1  Introduction 

In  this  chapter  we  will  develop  a  theory  that  will  enable  us  to 
propose  an  algorithm  for  solving  the  open-loop  optimal  dynamic  routing 
problem  for  single  destination  networks  with  unity  weightings  in  the  cost 
functional . 

In  Section  4.2  the  theoretical  background  is  developed,  and  in 
Section  4.3  the  algorithm  is  presented  (Section  4.3.1)  and  illustrated 
by  examples  (Section  4.3.2). 

Single  destination  networks  with  unity  weightings  in  the  cost  func¬ 
tional  have  properties  that  permit  the  solution  of  the  optimal  dynamic  open- 
loop  routing  problem  by  using  linear  programming. 

The  necessary  theory  for  developing  the  algorithm  will  be  based  on 
properties  from  Chapter  2,  also  the  cost  function  representation  in  the 
Z-plane  will  be  used  (see  Figure  3.1).  Since  it  is  known  that  there  always 
exists  an  optimal  solution  with  piecewise  constant  controls,  we  shall  search 
only  for  an  optimal  solution  with  this  property. 

Unity  weightings  in  the  cost  function  imply  a  *  1,  therefore 

i 


T*  p  • 

«  5  ■  l  xi 

i*l  1 


T 

«  X 


.T 

1  X 


n 

I* 

i*l 


(4.1) 


From  the  linear  program  (2.20)  and  from  Proposition  2.3  follows  that 
for  single  destination  networks  with  unity  weightings  in  the  cost  functional. 


the  optimal  solution  is  a  solution  of  the  following  linear  program: 


n 


subject  to 

and 


min  1  x. (t) 
i-1  x 


*i±0’  ^xieIp 


x.  =  0  ,  Vx.  e  B  , 
i  i  p 


then,  since  £  x. (t)  is  the  slope  of  the  solution  in  the  Z-plane  (see 


i=l 


Figure  4.1): 


Proposition  4.1 

The  slope  in  the  Z-plane  of  the  optimal  solution  is  minimal  among 
all  possible  slopes  satisfying 

x,  <  0  ,  V  x.  e  I 

l  -  '  l  p 

x.  =  0  ,  Vx.  e  B 
1  x  p 

Observe  that  since  the  slope  is  always  negative,  "minimal  slope"  means 
"the  steepest". 


Figure  4.1:  The  slope  in  the  Z-plane. 


□  Proposition  4.1 


Gpo.e^Hy  we  will  represent  trajectories  (x(t)  ,  te  )  in 

the  Z-plane,  since  this  allows  us  to  express  the  optimal  cost  as  the  area 
under  the  curve  representing  the  trajectory. 

Note  that,  for  a  given  time,  a  point  in  the  Z-plane  represents  all 

n 

the  state  vectors  with  the  same  value  of  J  x. . 

i=l  1 

We  conclude  from  Proposition  4.1  that: 


Conclusion  4.1 

In  single  destination  networks  with  unity  weightings  in  the  cost 
functional  there  is  no  feasible  trajectory  with  initial  slope  in  the  Z-plane 
smaller  than  the  slope  of  the  optimal  solution  (see  Figure  4.2). 


□  Conclusion  4.1 


4.2  Theoretical  Background 

Referring  to  a  feasible  solution  with  A  switches  (t  <  f  -  1  where 


f-  1  is  the  number  of  switches  in  the  optimal  solution),  we  define: 
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Definition  4.1 

The  feasible  solution  with  minimal  cost  function  among  all  the 
feasible  solutions  with  l  (or  less)  switches  will  be  called  the  "i-switches- 
opt ima 1 - so 1 ut ion" . 

□  Definition  .4.1 

The  algorithm  to  be  developed  in  Section  4.3  is  based  on  the  follow¬ 
ing  lemmas  and  theorems: 

Lemma  4 , 1 

If  there  exists  a  one-switch  trajectory  passing  from  xq  to  x^  in 
a  time  period  then  there  exists  a  non-switch  trajectory  connecting 

these  two  state  vectors  in  the  same  time  period  (see  Figure  4.3). 


Proofs 

We  will  denote  the  controls  in  the  one-switch  trajectory  as  follows: 

1  2 
u^ >  i  *  1 , . . .  ,m  for  t£  [^.tj),  and  u^  i«l,...,m  for  te[tlftf] 

(see  Figure  4.3). 
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The  trajectory  received  from  the  application  of  these  controls  is 
feasible,  therefore 

1  2 

u.  <  C.  and  u.  <  C.  ,  i  =  l,..,,m  .  (4.2) 

l  —  l  i—i 

We  now  refer  to  the  trajectory  received  from  the  application  of  the 


controls 

1  2 
U. T.  ♦  U.  T_ 

u  =  ,  i  =  l,...,m  ,  (4,3) 

x  Ti  +  t2 

during  the  time  period  t^,  (t^  =  Tj+Tj)  starting  from  the  initial  state 

vector  x  , 

-o 


Feasibility  Test 

(a)  >_  0,  Vi  =  l,..,,m  since  (4.3)  implies  that  is  a  sum  of  posi 
tive  elements. 

(b)  From  (4.3)  it  is  obvious  that 


u. 

l 


max{u*,u2}r,  +  max{u},u2>T_ 

_ l  i  1 _ l  l  2 

ri  +  T2 


then 


therefore, 


u.  <  max(u*,u2} 
i  —  i  1 


from  (4,2) 


% 


9 


u. 

1 


x  —  1,..., m  . 


Value  of  the  State  Vector  at  t. 


Since  x  is  a  linear  combination  of  the  elements  of  the  control 
vector  (see  (2.29)),  it  follows  from  (4.3)  that 


-1T1  *  -2T2 
T1  +  T2 


(4.4) 


•  * 

where  x^  and  x^  are  the  state  rates  of  the  one-switch  trajectory  before 


and  after  the  switch  respectively. 

At  t.  the  state  is  x  +  xt_;  then,  from  (4.4)  follows  that  the 
r  -o  -  i 

•  • 

state  is  x  +x,t  +  x.t  ,  i.e.  x.  (see  Figure  4.3). 

-o  -11  -l  £  -r 

In  light  of  these  facts, we  see  that  the  controls  defined  in  (4.3) 
determine  a  non-switch  trajectory  connecting  xq  and  x^  in  a  time  period  t 

□  Lemma  4,1 

Based  on  Lemma  4.1,  the  following  theorem  can  be  proved: 

Theorem  4.1 

If  there  exists  a  trajectory  with  f-1  switches,  taking  xq  to 
f 

xf  in  time  tf  (tf  =  \  -r),  then  there  exists  a  non-switch  trajectory 

connecting  these  two  state  vectors  in  the  same  time  period. 

Proof 

We  will  prove  the  theorem  by  repeated  use  of  Lemma  4,1. 

Suppose  that  there  are  f-1  switches  in  the  trajectory.  Then 
denote  by  x^x^, . , ,  ,x^  ^  the  state  vectors  at  tj,,..,t^  ^  respectively. 

Based  on  Lemma  4,1,  we  construct  a  non-switch  trajectory  passing 
from  xq  to  x2  in  a  time  period  T1  +  T2‘  With  this  tra0ector>r>  by 
using  again  Lemma  4.1,  we  construct  a  non-switch  trajectory  passing  from 
xq  to  x3  in  a  time  period  ti  +  t2  +  t3’  executing  this  procedure 
f  times,  we  get  a  non-switch  trajectory  passing  from  xq  to  x^  in  a 
time  period  t^  (see  Figure  4,4), 
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l 


□  Theorem  4 . 1 


We  conclude  from  Theorem  4.1  that: 

Conclusion  4_2 

If  for  a  given  initial  condition  x  there  exists  a  feasible  solution 

-o 

that  empties  the  network  (x^  =0)  in  a  time  period  t^,  then  there  exists  a 
feasible  solution  without  switches  that  empties  the  network  in  the  same  period, 

□  Conclusion  4.2 

All  remaining  theorems  in  this  chapter  refer  to  single  destination 
networks  with  unity  weightings  in  the  cost  functional. 


We  will  prove  first  the  following  fundamental  theorem: 


Theorem  4.2 


Any  feasible  trajectory  x(t)  must  satisfy 

n  n 


I  x  (t)  >  l  x*(t)  , 

i»l  1  i=l  1 


for  all  t  e  [tQ,t*]  • 


Proof 

We  can  put  in  words  the  statement  of  the  theorem  as  follows:  tl 
are  no  feasible  trajectories  in  the  Z-plane  penetrating  into  the  region 
underneath  the  optimal  trajectory.  Or  by  using  Theorem  4.1:  there  are  no 
feasible  trajectories  without  switches  in  the  Z-plane  penetrating  into  the 
region  underneath  the  optimal  trajectory. 

We  will  use  the  following  notation  (see  Figure  4.5): 
n 

(a)  X(t)  -  I  x.  (t) 

i=l  1 

n 

(b)  X  =  I  x. (t  ) 

o  1  o 

(c)  x^  =  (destination  node} 

We  will  prove  the  theorem  by  proving  that  there  is  no  trajectory 
without  switches,  x(t),  such  that  X(t')  <  X*(t')  for  some  t*e  [tp.tfl* 

We  shall  first  prove  the  theorem  for  the  case  when  the  optimal 


solution  has  at  most  one  switch. 
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n 


Lemma  4 . 2 

For  optimal  solutions  with  at  most  one  switch.  Theorem  4,2  holds. 

Proof 

From  Conclusion  4,1  follows  that  if  there  are  no  switches  in  the 
optimal  solution,  then  X(t)  >_X*(t),  te  since  the  initial 

slope  in  the  2-plane  of  the  optimal  solution  is  the  smallest  feasible  slope. 

Suppose  now  that  there  is  one  switch  in  the  optimal  solution;  from 
the  above  follows  that  X(t)  >_  X*(t)  for  te  t*) .  That  is,  if  there 
exists  a  time  t'  satisfying  X(t')  <  X*(t*),  then  t'  e  [t*,t*]  (see 
Figure  4.6).  Hence,  we  will  prove  the  lemrap  by  proving  that  X(t')  ^X*(t'), 

Vfe  [tJ,t*J. 


We  shall  refer  only  to  optimal  solutions  that  change  the  slope  in 

the  Z-plane  at  t*,  since  for  optimal  solutions  that  do  not  change  the 

slope  at  t*  a  new  optimal  solution  having  no  switches  can  be  constructed 

(by  direct  application  of  Theorem  4.1).  Note  that  a  solution  that  has  a 

switch  and  maintains  a  constant  slope  in  the  Z-plane  is  ;  solution  that 

n 

changes  the  controls  but  maintains  ][  x.  without  change. 

i=l  1 

From  the  above  follows  that  there  is  a  set  of  nodes  i  for  which 
(t)  vanish  at  t^  and  remains  zero  afterwards,  since  otherwise  it  would 
be  possible  to  continue  after  t*  with  the  same  slope  in  the  Z-plane, 
contradicting  the  assumption  that  there  is  a  switch  at  tj. 

We  denote: 

M  =  (x^/x?(t)  *  0,  Vte[t*,t*],  i  =  1, . . .  ,n>  , 

that  is,  M  includes  all  the  state  variables  in  the  optimal  solution  that 
are  zero  from  t  or  from  t*. 

O  1 

K  ■  {xi/x*^M,  i  ■  l,...,n)  , 
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that  is,  K  includes  all  the  state  variables  in  the  optimal  solution  that 
vanish  not  earlier  than  t^. 

Using  this  partition  of  the  network  we  define: 

Definition  4.2  (see  Figure  4.7) 

M  “Md")  ■  X  “jd'*)  ; 

(j  *d)eL 

UKd<«  ’  J,  \d">  - 

(k,d)eL 


(e)  U  (t)  =  l  I  u,  i  (t)  '  I  I  uiv(t)  ' 
keK  jeM  keK  jeM  JK 

(k ,  j )  eL  (j  »k)eL 


Observe  that  U^  is  the  net  flow  from  K  to  M,  i.e.  the  flow  from 
K  to  M  minus  the  flow  from  M  to  K. 

A 


(d)  C 


Kd 


l  c. 

keK 

(k,d)eL 


kd 


CWJ  and  CV1,  are  defined  in  a  similar  way: 

Md  KM 

(e)  XK(t)  *  l  x^t)  , 
keK 

XjdCt)  *  l  x  (t)  . 

"  icM  1 

•  •  • 

Likewise,  we  define  x(t),  X^Ct)  and  X^t)  as  the  time  derivatives  of 
X(t),  XK(t)  and  X^t)  respectively. 

<«  *0 • 

XCM  ^  ' 


□  Definition  4.2 
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We  will  prove  that  X(t')  >_  X* Ct ' ) ,  Vt*  e  [t* , by  developing 
appropriate  expressions  for  X(t')  and  X*(t*),  But  first  we  show  that 
X*  is  minimal  and  constant  on  [tQ,tp.  From  Figure  4.7  it  is  obvious 
that 

=  "u*Kd  -  um  *  (4-5) 

•  * 

hence,  in  order  to  show  that  XR  is  minimal  and  constant  we  will  show  that 

U*.  and  U*  are  maximal  and  constant. 

Kd  KM 

The  nodes  that  belong  to  K  vanish  only  at  t£,  therefore  the  flow 
of  information  from  these  nodes  to  the  destination  node  is  always  maximal, 
then 

UKd^  “  CKd  *  ^te[t0*tfJ  •  (4-6) 


For  each  node  that  belongs  to  M  the  net  flow  of  information  from 
K  to  M  from  t  until  the  switch  time,  is  smaller  than  or  equal  to  the 
flow  of  information  from  the  node  to  the  destination  node  (x  <_  0) .  There¬ 


fore,  from 


follows  that 


UKM^  <  UMd^  ’  te  [to'V 


(4.7) 


Likewise,  it  is  obvious  that  U*  (t)  for  t e  [t  ,t*)  is  maximal,  since 
otherwise  it  would  be  possible  to  increase  it  and  in  this  way  to  obtain  a 
trajectory  such  that  X(t)  *  X*(t)  for  te  [tQ and  X^t*)  >  °>  *■•••» 
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it  would  be  possible  to  obtain  a  new  trajectory  having  a  slope  smaller  than 
the  slope  of  the  optimal  solution  after  t*,  contradicting  Proposition  4.1. 

After  the  switch,  X^  is  equal  to  zero,  therefore  maximal  rate  of 

evacuation  of  the  information  from  the  network  means  maximal  X^.  Hence 

U*. C"t)  continues  to  be  maximal  after  the  switch.  Note  that  the  same  maxi- 
KM 

mal  flow  for  t  e  [tQ,t*)  is  possible  since  -  as  explained  above  -  this  flow 

was  feasible  when  x^  <_  0  for  i  e M,  then  it  is  obviously  feasible  when 

• 

it  is  needed  that  x.  =  0  for  ieM. 

l 

From  the  above  follows  that  U*  .(t)  is  constant  and  maximal  on 

KM 

[t  ,t^].  Therefore,  it  follows  from  (4.5)  and  (4.6)  that  X*  is  constant 

and  minimal  on  [t  , ttl . 

o  f 

We  will  use  now  this  fact  for  proving  the  lemma. 

From  the  definition  of  X  follows  that 


X=  XK+XM 


(4.8) 


likewise,  it  is  known  that  if  x  is  constant  on  [tQ,t'] 


x(t’)  =  x  +■  xt’  , 
-o 


(4.9) 


where  without  loss  of  generality  it  is  assumed  that  t  *  0. 

o 

•  4th  *  ^ 

Using  the  facts  that  =  0  on  [t^tj  and  XR  is  constant  on 
[0,tp,  it  follows  from  (4.8)  and  (4.9)  that 


where 


xoK+*;t,+xoM+ifc;  * 


X  *  X  u  +  X  „ 
o  oM  oK 


(4.10) 


But  the  contents  of  nodes  that  belong  to  M  vanish  at  tj,  then 


^  -XoM  +  XK  =  0 


(4.11) 


By  using  (4.11),  equation  (4.10)  may  be  written  as 

X*(f)  -X0K*^t'  ,  t'MtJ.t*]  •  C4.12) 

The  trajectory  we  are  going  to  test  has  no  switches,  therefore 

■  XoK*V'*XoM*V'  *  C4  l3) 

but 

Vt,J  ■  XoM*V  '  t4’14) 

then,  (4.13)  may  be  written  as 

X(t')  =»  XoK  +  XKt' ♦XM(t')  ,  fe[tj,t*].  (4.15) 

•  •  •  a 

0  and,  since  X*  is  minimal,  X*  <_  X^.  There¬ 
fore  comparing  (4.15)  with  (4.12)  we  conclude  that 

X(t')  >  X*(f)  ,  t»e[tj,tj]  . 

□  Lemma  4 . 2 

We  will  use  now  Lemma  4.2  in  order  to  prove  by  induction  the  theo¬ 
rem  for  the  general  case.  For  this  purpose  we  will  suppose  that  Theorem  4.2 
holds  for  an  optimal  solution  having  h-  1  switches,  and  then  we  prove  that 
the  theorem  also  holds  for  an  optimal  solution  having  h  switches. 

For  an  optimal  solution  with  h  switches,  there  cannot  be  penetra¬ 
tion  in  the  Z-plane  into  the  zone  underneath  the  optimal  trajectory  for  any 
time  smaller  than  t*  (because  the  theorem  holds  for  an  optimal  solution 

a 

having  h-1  switches),  i.e.,  X(t)  ^  X*(t)  for  te[tQ,t*)  where  t^ 

A 

is  the  h-switch  in  the  optimal  solution.  It  follows  then  that  X (t )  can 
be  smaller  than  X*(t)  only  on  [t?,t£]  (see  Figure  4.8). 


It  is  known  that  X^t’)  > 


Figure  4.8:  Trajectories  in  the  Z-plane  for  proving 
Theorem  4.2  in  the  general  case. 

The  proof  will  follow  the  same  pattern  as  the  proof  of  the  lemma. 
Therefore  we  define  the  sets  K  and  M,  then  we  show  that  X*  is  minimal 
and  constant,  and  as  a  consequence  of  this,  that  X(t')  is  greater  than  or 
equal  to  X*(t'). 

We  then  define: 

K  =  {set  of  nodes  whose  contents  vanish  in  the  optimal  solution 
exactly  at  t*}  , 

M  =  {set  of  nodes  whose  contents  vanish  in  the  optimal  solution 
before  t^}  , 

that  is,  the  set  M  includes  all  nodes  that  are  empty  from  tQ,  and  those 
that  are  empty  from  t*,t*, .. .,t£. 

In  order  to  prove  that  X*  is  minimal  and  constant,  suppose  that 

*  * 

XK  is  not  minimal  in  some  time  interval.  It  is  then  possible  to  increase 
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the  net  flow  from  K  to  M  in  this  interval.  As  a  consequence  of  that,  a 
trajectory  that  coincides  with  the  optimal  trajectory  in  the  Z-plane  until 

^)  >  0  is  obtained.  With  this  new  trajectory  it  is  pos¬ 
sible  to  continue  after  t^  with  a  slope  smaller  than  that  of  the  optimal 
trajectory,  contradicting  Proposition  4.1.  X*  is  constant  in  [t  ,t*] 
because  of  the  fact  that  the  same  maximal  flow  from  K  to  M  remains  before 
and  after  every  switch  (note  that  at  each  switch  some  of  the  state  variables 
that  belong  to  M  become  zero) .  This  is  explained  by  the  fact  that  if  this 
maximal  flow  was  possible  before  these  states  become  zero  (i.e.,  when  in 
those  nodes  the  derivative  of  the  state  was  negative),  it  is  obvious  that  it 
is  still  possible  after  the  switch  when  these  nodes  are  empty  (i.e.,  when 
the  derivative  of  the  state  is  equal  to  zero) . 

By  making  the  same  steps  as  in  Lemma  4.2,  the  same  expressions  for 
X*(t')  and  X(t')  are  obtained,  that  is: 

x*(t’)  =xoK  +  x*t. 

and 

*(»’>  *  XoK*V*Vt')  • 

A  •  A  £ 

Since  X^t*)  is  nonnegative  and  XR  ^  X^,  the  theorem  is  proved  by  compar¬ 
ing  these  two  expressions. 

□  Theorem  4 . 2 

Based  on  Theorem  4.2  the  following  theorems  can  be  proved: 

Theorem  4 , 3 

Referring  to  solutions  with  piecewise  constant  controls.  Each  optimal 
solution  to  the  minimum  delay  problem  is  also  optimal  to  the  minimum  total 
time  problem  (i.e.,  min  t^) . 


t£  such  that  X^Ct 
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Proof 

From  Theorem  4.1  follows  that  if  there  exists  a  trajectory  with 
total  time  smaller  than  t^  ,  then  there  exists  also  a  non-switch  trajectory 
with  total  time  r  smaller  than  t^  ,  i.e., 

t  <  tj  .  (4.16) 

Proposition  4.1  implies  that  the  trajectory  satisfying  (4.16)  has  an  initial 

slope  greater  than  the  slope  of  the  optimal  solution.  Therefore,  the  ful- 

n  n 

fillment  of  (4.16)  requires  that  £  x.(t')  <  £  xf(t')  for  some  t' 

i=l  1  i=l  1 

(where  x(t)  is  the  trajectory  satisfying  (4.16)),  contradicting  Theorem  4.2 
(see  Figure  4.9).  Hence,  t^  is  the  optimum  of  the  minimum  total  time 
problem. 


X 


,  □  Theorem  4,3 

Recall  that  we  deal  with  single-destination  networks  with  unity 
weightings  in  the  cost  functional. 


Theorem  4.4 


A 

Let  x(t)  be  the  one-switch-optimal -solution  and  t*  its  switch 
time.  Then  we  have 

X(t’)  =  X*  (t ' )  . 

Or  in  words,  the  one-switch-optimal-solution  represented  in  the  Z-plane  has 
its  switch  exactly  on  the  optimal  trajectory. 

Proof 

We  will  prove  the  theorem  by  constructing,  for  any  one-switch  solu¬ 
tion  whose  switch  point  is  not  on  the  optimal  trajectory  in  the  Z-plane,  a 
one-switch  solution  with  a  switch  point  exactly  on  the  optimal  trajectory 
and  with  a  smaller  delay  (see  Figure  4.10). 

Consider  an  arbitrary  one-switch  trajectory  x(tj  whose  switch  point 
is  not  on  the  optimal  trajectory  in  the  Z-plane.  Then,  from  Theorem  4.2  follows 
that 

X(t')  >  X*(t') 

A  * 

and  from  Theorem  4.3  follows  that  t^  >_  tf  . 
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We  construct  now  a  one-switch  trajectory  x* (t) ,  whose  switch  point  is  on 
the  optimal  trajectory  in  the  Z-plane,  having  a  cost  smaller  than  the  cost 
of  the  one-switch  solution  considered  above.  Theorem  4.1  implies  that  there 
exists  a  non-switch  trajectory  taking  xq  to  x*(tr)  in  a  time  period  t', 
and  a  trajectory  taking  x*(t')  to  zero  in  a  time  period  t£-t*.  Let  x'(t) 
be  this  overall  trajectory.  This  constructed  one-switch  trajectory  has  a 
cost  smaller  than  the  cost  of  x(t)  since  they  have  the  same  switch  time, 
X'(t')  <  X(t')  and  t£  =  t*  <_  tf  (see  Figure  4.10). 

□  Theorem  4 . 4 


Definition  4.3 


is  the  average  slope  in  the  Z -plane. 


Definition  4.4 


X. 

l 


aXCV-xcvj) 

Wi 


that  is,  )L  is  the  slope  of  a  trajectory  represented  in  the  Z-plane  during 


"i-i’V- 


Theorem  4.5 

The  switch  time  in  one -switch-optimal -trajectories  occurs  at: 

(a)  a  switch  time  t?  of  the  optimal  trajectory,  if  X^  <  D  and 

>  D> 

(b)  or  at  any  time  in  the  time  interval  [t±_1 , ,  if  Xi  =  D, 

*I+i  > D  and  K-i <  D* 


Case  (a) 


v. 


Case  (b) 


Figure  4.11:  The  two  cases  of  Theorem  4.5 
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PTOOf 

We  will  find  the  switch  time  t’  of  the  one-switch-optimal -trajectory 
by  representing  the  cost  function  J  as  a  sole  function  of  t'  and  then, 
by  differentiation,  we  will  find  the  optimal  t'. 

The  cost  function  for  a  one-switch  trajectory  will  be  represented 
as  the  sum  of  the  areas  of  two  triangles  (see  Figure  4.12). 


X 


From  Figure  4.12  follows  that 


J  = 


X(t')t. 


For  the  one-switch-optimal-trajectory  t^  ■  t^  .  f°T  this  case 


X  t'  X(t')t* 
o  ^  f 


(4.17) 


From  (4.17)  we  see  that  J  is  determined  by  t'  and  X(t').  But 
in  Theorem  4.4  was  proved  that  X(t')  is  located  exactly  on  the  optimal 
trajectory  in  the  Z-plane,  therefore  X(t')  can  be  written  as 


X(t’)  =  X0  +  Vl  +  X2T2  +  •••  +  VlTl-l+Xi(t,"V  » 


where  t  -  1  is  the  number  of  switches  in  the  optimal  solution  from  t 
to  t',  and  Xj  xs  the  sum  of  the  elements  of  the  state  rate  vector  in 


[t?  .,t?).  Then 

j-1  y 


J(t’) 


xt-  v 

0  j*l  J  3 


t*{’  t^<f<t;.  (4. 


By  differentiating  (4.19)  with  respect  to  t'  we  get 


X  X*t* 
-2.  +  1  £ 
2  2 


and  by  replacing  D  (see  Definition  4.3)  into  (4.20)  we  have: 


★ 

dJ  -  tf 
dtr-  (-D  + V  T 


For  continuing  the  proof  we  need  the  following  proposition: 


Proposition  4.2 


The  optimal  solution  satisfies 


*;  i  *;+1  >  3  =  1 . f-i 


Proof 


From  Proposition  2.3  follows  that  the  optimal  solution  satisfies 
that  X(t)  is  minimal  Vt  subject  to  the  constraints  xi  <^0,  V  t  1^ 
and  x^  *  0,  Vx^  e  B^.  Since  we  deal  with  optimal  solutions  with  piece- 
wise  constant  controls  satisfying  x  <_  0,  the  number  of  empty  nodes  can 
only  be  increased  at  switching  times.  Therefore  some  of  the  constraints 
of  the  type  <_  0  can  be  transformed  after  a  switching  time  to  con¬ 
straints  of  the  type  x.  ■  0. 
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Hence,  after  a  switching  time,  the  constraint  region  can  be  equal  to 
or  smaller  than  the  constraint  region  before  the  switching  time,  and  hence 


Xj  ~  X3+l’  j 


□  Proposition  4.2 


From  Proposition  4.2  follows  that  if  X*  <  D,  then  X?  <  D,  Vj  <  Jt, 


and  hence,  from  (4.21),  gp  <0,  V t •  e  [t  ,t^] . 


Now  let  t*  be  the  switching  time  of  the  optimal  solution  for  which 
• ,  •* 

X^  <  D  and  Xi+^  >  D.  Then  it  follows  that  the  switching  time  t*  of  the 
one-switch-optimal -solution  is  t?  (see  Figure  4.11,  case  (a)). 

On  the  other  hand,  it  is  possible  that  there  exists  a  time  interval 
[t?  ^,t?)  in  which  X?  =  D;  then  from  Proposition  4.2  follows  that 
X?  j  <  D  and  X*+1  >  D,  and  the  switching  time  t'  of  the  one-switch- 
optimal-solution  satisfies  t' e  [t?  ^,t?)  (see  Figure  4.11,  case  (b) . 


□  Theorem  4 . 5 


Theorem  4.6 

—»■■■■■  ■■■  ■  • 

If  the  one-switch-optimal-solution  has  no  switches,  then  the 
t-switches-optimal-solution  (i>l)  has  no  switches. 

Proof 

We  will  prove  the  theorem  by  contradiction. 

Suppose  that  there  exists  an  t-switches-optimal  trajectory  with 
a  cost  smaller  than  the  cost  of  the  non-switch  solution.  From  Theorem  4.3 
follows  that  in  these  two  solutions  the  final  time  is  t£  (since,  if  the 
final  time  is  greater  than  t£,  then  it  is  possible  to  decrease  the  cost 
contradicting  the  optimality  assumption) . 


The  above  fact  and  Proposition  4.2  imply  that  all  switches  of  the 
i-switches-optimal -trajectory  are  in  the  Z-plane  underneath  the  non-switch- 
optimal  -trajectory.  Then,  Theorem  4.1  implies  that  it  is  possible  to  con- 
struct  a  one-switch -optimal -trajectory  with  a  smaller  cost  than  the  non- 
switch-trajectory,  contradicting  the  assumption  that  the  one-switch-optimal 
solution  has  no  switches. 


X 


4.3  Algorithm  for  Solving  the  Open-Loop  Optimal  Dynamic  Routing  Problem 
for  Single-Destination  Networks  with  Unity  Weightings 
in  the  Cost  Functional 

4.3.1  The  Algorithm 

We  present  now  an  algorithm  that  solves  the  open-loop  optimal  dynamic 
routing  problem  for  single-destination  networks  with  unity  weightings  in  the 
cost  functional. 


The  proposed  algorithm  is  based  on  the  theory  developed  in  Section  4.2 
We  present  now  a  brief  description  of  the  algorithm,  and  then  we 


explain  it  in  detail.  At  the  end  of  this  section  the  algorithm  is  presented 
as  composed  of  two  operations,  with  the  second  one  applied  repeatedly. 


-  61  - 


The  algorithm  finds  an  optimal  solution  by  solving  a  series  of  linear 
programs.  This  is  done  by  finding  first  the  non-switch-optimal-solution  and 
the  final  time  of  the  optimal  solution,  i.e.  t*.  Then  a  switch  is  allowed 
in  the  solution,  so  that  the  algorithm  finds  now  the  one-switch-optimal - 
solution.  The  knowledge  of  t£  allows  to  calculate  the  one-switch-optimal- 
solution  by  solving  a  linear  program.  The  algorithm  progresses  by  finding 
the  3,  7,  15,  etc. -switches -optimal -solutions  by  allowing  one  extra-switch 
in  each  of  the  intervals  of  the  previous  solution  until  the  optimal  solution 
is  found.  In  each  step,  the  solutions  are  obtained  by  solving  a  linear 
program . 


Explanation  of  the  Algorithm 

From  Theorem  4.3  follows  that  the  final  time  of  the  optimal  solution 
is  t*,  where  t*  is  the  minimal  final  time  among  all  feasible  solutions 
having  piecewise  constant  controls.  This  time  is  calculated  by  using  Theo¬ 
rem  4.1  which  says  that  for  any  trajectory  passing  from  xq  to  x^  in  a 
time  period  t^,  there  exists  a  non-switch  trajectory  passing  from  xq  to  x^ 
in  the  same  time  period.  Hence  t£  will  be  found  by  looking  for  the  non¬ 
switch  trajectory  having  minimal  final  time,  that  is  by  solving  the  linear 
program: 


t^  =  min  t 

subject  to 

x  +  Bz  =  0 
-o 

Dz  -  Ct  <_  0 
z  >  0,  t  >_  0 


(4.22a) 


(4.22b) 


Note  that  the  constraints  (4.22b)  are  the  constraints  (3.24)  -  (3.27)  for  a 
non-switch  trajectory. 


The  knowledge  of  t^  allows  us  to  obtain  the  one-switch-optimal- 
solution  by  solving  a  linear  program.  In  order  to  do  so  we  represent  the 
cost  function,  for  a  one-switch  trajectory,  as  the  sum  of  the  areas  of  two 
triangles  (see  Figure  4.12  and  equation  (4.17)). 

T 

If  we  denote  by  t^  the  switch  time,  then  X(t^)  =  -1  Bz^  (see 
Figure  4.14) . 

X 


t 


Figure  4.14:  The  performance  index  J  for  the  one-switch-optimal-solution 
Hence,  from  (4.17)  follows  that 


i  -2  f 

J  ’  2  '  2 

Therefore,  the  one-switch-optimal  solution  can  be  obtained  by  solving  the 
linear  program: 
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XT,  1TBz 

min  J  »  t *  ,  (4.23a) 

subject  to 

x  +  Bz  >  0 
-o  -1  — 

x  +  Bz,  +  Bz. 

-0  -1  -L 

Dz.  -  Ct .  <  0 

-l  -  i  — 

z.  >  0 

-i  — 

T.  >  0 
l  — 

Note  that  the  constraint  region  (4.23b)  is  the  constraint  region  (3.24)  - 
(3.27)  for  the  one-switch  case. 

If  the  optimal  solution  has  no  switches,  there  will  be  no  decrease 
in  the  cost  function  in  this  step  of  the  algorithm.  If  the  optimal  solution 
has  exactly  one  switch,  then  it  is  received  in  this  step  of  the  algorithm 
and  in  the  next  step  there  will  be  no  decrease  in  the  cost  function  (this 
fact  will  be  explained  later). 

From  Theorem  4.4  follows  that  X'(tp  is  equal  to  X*(tp  where 

x'(t)  is  the  one-switch-optimal-trajectory  and  tj  its  switch  time.  Note 

that  X'(tp  is  indeed  equal  to  X*(tp,  but  it  is  possible  for  x'  (tp 

to  be  different  from  x*(tp.  If  these  two  vectors  were  equal  it  would  have 

been  possible  to  separate  the  problems  into  two  problems:  a  problem  with 

initial  condition  x  ,  end  condition  x'(t')  and  final  time  t ' ,  and  a 
-o  -  1  i 

problem  with  initial  condition  x'(tp,  end  condition  zero  and  final  time  t£. 
Since  this  is  not  the  case  however,  the  problem  is  more  complicated.' 

With  the  results  of  the  results  of  the  one-switch-optimal -solution 
as  a  basis,  we  will  find  next  the  best  solution  among  all  trajectories  pas¬ 
sing  through  X'(tp  and  having  two  more  switches  in  addition  to  the  switch 


•  •* 
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at  tj,  a  switch  in  and  another  in  [t^,t^]..  This  problem  can  be 
solved  using  linear  programming  if  we  represent  the  cost  function  as  the  sum 
of  the  areas  appearing  in  Figure  4.15. 


X 


Figure  4.15:  Performance  index  as  the  sum  of  areas 
for  a  three-switches  solution. 


The  expression  obtained  is: 

tx0-x'(tpj 


+  x'ttpxj  + 


X'(tpt3 


1V2 


(4.24) 


where  the  number  under  each  expression  denotes  the  corresponding  area  in 
Figure  4.15.  All  the  variables  refer  to  the  three -switches  trajectory  and 
all  known  values  (from  the  previous  step  of  the  algorithm)  are  marked  with 
a  "prime". 
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The  constraints  (3.24)- (3.27)  will  be  written  in  two  parts,  'a  part 
before  t|  and  the  other  from  tj  until  t^.  In  addition  it  is  known  that 
at  t’  the  state  vector  x(tp,  which  must  be  determined,  satisfies 

!  x.(t')  =»  X'(t')  , 

i=l  1  1  1 

and  that  the  additional  constraints  T1  +  T2  =  anc*  T3  +  T4  *  T2  must 

* 

hold.  It  follows  then  that  the  constraint  region  is: 


xQ  +  Bz^  >_  x(tp 


xq  +  Bzx  +  Bz2  =  x(tp 


I  x  (tp  =  X'(t') 
i»l  1  1  1 


x(tp  +  Bz3  0 
x(tp  +  Bz3  +  Bz4  =  0 


T  +  T  =  T  1 
1  2  1 


T3  +  T4  =  T2 


Dz.  -  Ct.  <  0 

-l  -  l  — 


z.  >  0 
-l  — 


i  =  1, ...  ,4 


(4.25) 


T.  >  0  I 
1  —  ' 

Hence,  the  linear  program  to  be  solved  is  to  minimize  (4.24)  subject  to 


(4.25). 

If  there  is  a  decrease  in  the  cost  function,  a  switch  is  added  in 
each  of  the  time  intervals,  a  new  linear  program  is  solved,  and  the  procedure 
is  repeated  until  the  optimal  solution  is  found. 

The  algorithm  can  be  described  in  a  more  concise  way  if  we  look  at 


it  as  composed  by  two  operations,  as  follows: 


Find  the  final  tine  of  the  optimal  solution  (t^)  and  the  non-switch 
optimal -solution  by  solving  the  linear  program 


subject  to 


tj  ■  min  t 


x  ♦  Bz  ■  0 
-o 


Dz-Ct-0  (4.26) 

z  >  0 
t  >0 

Operation  2 

A  feasible  solution  x'(t)  with  £  switches,  satisfying 


X*(t«)  *  X*(t') 

for  all  its  switching  times  t?,  i  *  1,...,&,  is  given.  Find  the  feasible 
solution  with  minimum  cost  function  among  all  solutions  having  a  switch  in 
each  of  the  time  intervals  (tJ,t£+p,  i  *  0 (t£+1  “  tP*  in  ad(*i- 
tion  to  the  switches  t!,  i  *  1,...,£,  and  satisfying  X'(t!)  equal  to 
X*(tp,  i  =  1,...,A.  This  solution  will  have  21  +  1  switches  and  is 
obtained  by  solving  the  linear  program  (see  Figure  4.16): 


min  J 


1*0 


X»(tp  -X'(t'+1) 
- 5 - 


2i+l 


>  + 


% 

l 

i*l 


X'(t')x'  , 


(4.27a) 


subject  to: 


?0  +  B-1  - 

XQ  +  BZj  +  Bz2  -  X(tp 


y  x.(tn  -  x'ttn 

i»l  1  1  1 

x(tp  +  Bz3  >  x(tp 


x(t')  +  Bz3  +  Bz4  =  x(tp 


I  X.(t')  =  X'(t') 
i=l  1  1  1 


*  b;2i-i  *  bj2i  * 

*"i>  *  Bht.l  - 0 


(4.27b) 


?";>  *  Bht*i  *  b;2**2  *  ° 


T,  +  T  =  T  ' 
1  2  1 


T3  +  T4  =  T2 


T  +  T  a  T  * 

21+1  2)1+2  Vl 


Dz^  -  Ct^  0 


Zi  >_  0  i  =  1, . . .  ,2)1+2 


t.  >  0 

l  — 


Note  that  \  X’(t!)T.'  is  constant  (is  received  from  the  previous  step  of  the 
i*l  1  1 

algorithm)  and  appears  in  J  only  for  i.  >  1. 
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These  two  operations  are  applied  as  follows: 


Figure  4.17:  The  two  operations  of  the  algorithm. 

Next  we  shall  prove  that  the  algorithm  converges  to  the  optimum  in 
a  finite  number  of  steps.  In  order  to  do  so  we  shall  prove  first  the  follow¬ 
ing  propositions: 

Proposition  4.3 

Given  that: 

X'(tj)  -  X*(tj),  j  =  1.....4  , 

where  x’ (t)  is  the  solution  of  Operation  2  with  l  switches.  Then 

X(t.)  =  X*(t.),  i  =  1,...,2*+1  , 


where  x(t)  is  the  solution  of  the  next  step  of  the  algorithm,  i.e.  x(t)  is 
the  solution  of  Operation  2  with  2£+l'  switches. 
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Proof 

From  Theorem  4.2  follows  that  X(t^)  >_X*0~),  i  ■  l,...,2g+l. 

Hence,  if  the  proposition  is  not  satisfied,  X(t^)  >  X*^)  for  at  least 
one  switching  time. 

We  now  assume  that  there  exists  a  solution  to  Operation  2  which  does 
not  satisfy  the  proposition,  and  then  construct  a  solution  to  Operation  2 
that  does  satisfy  the  proposition,  having  a  smaller  cost  function. 

The  constructed  solution  (see  Figure  4.18)  is  the  solution  passing 
through  x*(t^)  for  all  i  =  l,...,2i+l.  It  is  clear  that  for  this  solution 
X(t^)  is  equal  to  X*(t^)  at  all  switching  times,  and  its  cost  is  smaller 
than  the  cost  of  the  assumed  solution. 


X 


Figure  4.18:  Trajectories  in  the  Z-plane  for  Proposition  4.3 


□  Proposition  4 . 3 


71 


Proposition  4,4 


If  after  a  step  of  the  algorithm  there  is  no  decrease  in  the  cost 
function,  then  further  steps  (allowing  more  switches  in  each  time  interval) 
will  not  decrease  the  cost  function. 


Proof 


A  step  of  the  algorithm  means  allowing  a  switch  in  each  of  the  time 


intervals  [t!,t|+1],  i  *  0 where  tJ,  i  =  1,..,,A  are  the  switch¬ 
ing  times  obtained  in  the  previous  step  of  the  algorithm.  Suppose  then, 
that  the  present  step  does  not  decrease  the  cost  function  value,  and  look  at 
the  problems  defined  by  the  switching  times,  that  is,  look  at  the  problems 
with  initial  condition  x(tp,  end  condition  x(t^+1),  and  final  time  rj, 
i  =  l,..,,£fl.  From  Theorem  4.6  follows  that  allowing  more  switches  in  each 
of  these  problems  will  not  decrease  the  cost  function.  Hence,  further  steps 


will  not  decrease  the  cost  function. 


Theorem  4,7:  Convergence  Theorem 


□  Proposition  4,4 


The  proposed  algorithm  converges  to  the  optimal  solution  in  a  finite 


number  of  steps. 


Proof 


The  optimal  solution  has  a  finite  number  of  switches  distributed  in 
[tQ,t*] .  In  each  step  of  the  algorithm  the  allowed  number  of  switches  in¬ 
creases,  in  such  a  way  that  in  a  finite  number  of  steps  we  will  obtain  a 
solution  with  the  required  number  of  switches  distributed  in  [to,tp,  like 
the  optimal  solution.  From  Proposition  4,3  follows  that  the  solution  of 
Operation  2  for  this  case  satisfies  X(t^)  .equal  to  X*(t^),  and  hence  it 
is  the  optimal  solution. 


□  Theorem  4 . 7 


The  theorem  above  implies  that  after  a  finite  number  of  steps  an 
optimal  solution  will  be  achieved.  From  Proposition  4.4  follows  that  the 
stopping  criterion  of  the  algorithm  is  no  decrease  of  the  cost  function 
value  after  some  step;  that  is,  if  in  a  certain  step  there  is  no  decrease 
in  the  cost  function  value,  the  optimal  solution  was  reached  in  the  pre¬ 
vious  step. 

The  algorithm  proposed  in  this  section  was  programmed  (see 
Appendix  C) ,  and  several  networks  were  solved.  The  program  has  the  property 
that  it  discovers  when  the  time  between  switches  is  zero,  and  for  this  case 
it  does  not  allow  more  switches  in  this  interval.  In  this  way  unnecessary 
new  variables  are  defined  in  the  following  steps  of  the  algorithm. 

In  the  next  section  we  present  some  of  the  networks  that  were  solved. 


4.3.2  Examples 

In  this  section  we  present  numerical  results  obtained  for  the  solu¬ 
tion  of  the  optimal  open- loop  routing  in  a  number  of  networks.  The  prob¬ 
lems  were  solved  using  a  computation  program  written  in  FORTRAN  (see 
Appendix  Cl.  The  program  was  run  on  an  IBM  370/168  computer  and  for  each 
example,  besides  Example  4,3,  a  maximal  memory  region  of  51 2K  was  used. 

We  present  now  three  examples .  In  the  first  we  present  the  network 
parameters,  the  results  of  the  different  steps  of  the  algorithm  and  the 
optimal  solution.  In  the  other  two,  only  the  network  parameters  and  the 
optimal  solution  are  printed. 

Example  4.1 


Refer  to  the  network  appearing  m  Figure  4.19  with  initial  condition 
T 

x,  -  (2,5,4)  ..  The  numbers  on  the  links  are  their  capacities. 
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The  results  of  the  program  are: 

THE  CAPACITY  VECTOR  ISJ 

2.000  1.000  2.000  1.000  1.000  1.000  4.000 


THE  B  MATRIX  IS: 

-1.0-1 .0-1.0  0.0  0.0  0.0  0.0 
1.0  0.0  0.0— 1.0— 1.0  1.0  0.0 
0.0  1.0  0.0  1.0  0.0-1.0-1  .0 

THE  INITIAL  CONDITION  ISi 
2.00  5.00  A. 00 


S  OPTIMUM  COST  VALUE  2.49999HE  OO 


t  ERROR  FLAGS 

0 

ANO  2.980232E-07 

:  OPTIMUM  SOLUTION 

1 

0.0 

2 

0.0 

3 

1.999998F  90 

4 

2.499996E  00 

5 

£‘.5u000l>E 

00  6 

0.0 

7 

6.4999'JOE  00 

8 

2.499900C 

00 

E 

TOTAL  TIME  I S! 0.24999M9F  01 

I  OPTIMUM  COST 

VALUE  I.I/2E000E 

Ol 

1  ERROR  FLAGS 

-1 

ANO  4.758372E 

-06 

:  optimum  solution 

1 

0.0 

2 

0.0 

3 

2.000000E 

00 

4 

1  .333333F  00 

5 

1.33J33/E  00 

6 

0.0 

7 

S.333332E  00 

b 

u.u 

9 

0.0 

10 

0.0 

11 

1.1 6666 YE  OO 

12 

1.1 6666 YE 

00 

13 

0.0 

14 

1.  16666FE  00 

1  S 

1  .  J333  3.*E 

00 

16 

1  . 1 66667F  00 

17  - 

l.I‘»209a:-06 

:  optimum  cost  valut  p.hossspe  oo 


S  ERROR  FLAGS  -1 

:  OPTIMUM  SOLUTION 

ANO  1.647744? 

-06 

1 

O.b 

2 

0.0 

3 

2.000000E  00 

4 

1.0000 OOE  00 

5 

1 .00000 OE 

00 

6 

0.0 

7 

4.000000E  00 

H 

0.0 

A 

0.0 

10 

0.0 

1  1 

3.333  336E-01 

1? 

3. 333335? -01 

13 

0.0 

14 

1.33333  IE 

00 

15 

0.0 

16 

0.0 

17 

0.0 

1H 

5.  <>6046 At —07 

19 

5 .960  4  64  E— 07 

20 

0.0 

21 

2.S4 3130E-06 

22 

0.0 

23 

0.0 

24 

0.0 

25 

1. 16666 7E  00 

26 

1. 166667C 

00 

27 

0.0 

28 

1.166667E  00 

29 

0.0 

30 

2  •  333333F  00 

31 

1 .0530 ISE— 06 

32 

l.OOOOOOE 

00 

33 

3. 33333VE— 01 

34 

37 

5.960464E-07 
8.145968!: -07 

35 

1. 166667E 

00 

36 

0.0 

:  OPTIMUM  COST  VALUE  5.1 3885 7E  00 


S  ERROR  FLAGS 

-1 

AND  2 .574  921 E 

-05 

1 

S  OPTIMUM  SOLUTION 
0.0  2 

0.0 

3 

1.999998E  00 

4 

9.999996E— 01 

5 

9.99999OE-0I 

6 

0.0 

7 

3.999990E  00 

e 

0.0 

9 

0.0 

10 

0.0 

ti 

9.53674  3E-07 

12 

9.536743E-07 

13 

0.0  __ 

14 

6.675 720E-06 

1  5 

0.0 

16 

0.0 

17 

1. 788 1 39E—06 

18 

3.333327E— 0 1 

19 

3.333327E-01 

20 

0.0 

21 

1 . 333328E  00 

22 

0*0 

23 

0.0 

24 

0.0 

25 

7.9A7273E-08 

26 

7.947273E-00 

27 

0.0 

28 

7.947273F-08 

29 

0.0 

30 

0.0 

31 

0.0 

32 

5.960464E-O7 

33 

S.960464E-07 

34 

0.0 

35 

o.o 

36 

0.0 

37 

0  .0 

3B 

0.0 

39 

0.0 

40 

1  •  9  07  348  E  *“06 

41 

0.0 

4  2 

0.0 

43 

0.0 

44 

0.0 

45 

0.0 

46 

9.1 39381 E— 07 

47 

9.  13938  IE -07 

48 

0.0 

49 

2.1060  30F.-06 

50 

0.0 

51 

0.0 

52 

0.0 

53 

1.1 666646  00 

54 

1.166664E  00 

55 

0.0 

56 

1.166664E  00 

57 

1.768139E-06 

58 

2 .999090*  00 

59 

9.040QBt£  -0 1 

60 

0.0 

61 

2.333333ft  00 

62 

I.152356E— 06 

63 

0.0 

64 

2.333330E  00 

65 

1.746403E-06 

66 

9  .99999  6E  —0 1 

67 

9.53674JE-07 

68 

3.  33332  7E -01 

69 

7.947273E-08 

70 

5  .960*  64E— 07 

71 

0.0 

72 

9.1 J9381E-07 

73 

1.166664E  00 

74 

0.0 

75 

0.0 

76 

0.0 

77 

2.264976E-06 

TAU (  1  ) =0  > 1 OOOOOOE  01 


U  <  1  )  =0  -0 

UC  21*0. 0 


ue  3 > *0 .200000 CE  01 

u<  *  1=0 .1 OOOOOOE  01 

U(  S ) *0 • 1 OOOOOOE  01 


U<  6  )  *0 .0 


U(  7 ) *0 .4000000E  01 


76 


TAOI  7 > *0.333 333 5E  00 


uc 

11*0.0 

uc 

?>*0.0 

UC 

3 ) *0  *0 

uc 

♦  >=0.1 OOOOOOE 

01 

U( 

f>)  =0.1  OOOOOOE 

01 

u< 

61=0.0 

U{ 

7)  =0.3900P«»ie 

01 

TAUC  31=0.1 166667E 

UC 

<1 

o 

* 

o 

UC 

71=0.0 

U( 

31=0.0 

uc 

4 1=0.1 OOOOOOE 

01 

uc 

5 ) =0 • 1 OOOOOOE 

01 

uc 

61*0.0 

UC  7) *0 • ! OOOOOOE  01 
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In  the  first  part  we  print  the  network  parameters,  i.e.  the  capacity 
vector,  the  B-matrix  and  the  initial  conditions.  The  first  element  of  the 
-capacity  vector  is  the  capacity  of  the  link  that  corresponds  to  the  first 
column  of  the  B-matrix,  and  so  on.  In  the  second  part  the  results  of  the 
different  steps  of  the  algorithm  are  presented  (the  meaning  of  the  different 
variables  of  the  linear  programs  appears  in  Appendix  C) .  Note  that  in  this 
example,  there  is  no  elimination  of  variables  when  a  time  interval  is  zero 
in  order  to  get  an  easy  reference  to  the  meaning  of  the  variables.  The 
steps  of  the  algorithm  are: 

(a)  Operation  1 :  In  the  first  linear  program,  the  final  time  (TOTAL 
TIME)  is  calculated,  and  the  non-switch-optimal-solution  is  obtained. 

(b)  Operation  2  for  1=0:  The  one-switch-optimal-solution  is  calcul¬ 
ated.  Since  there  is  an  improvement  in  the  cost,  this  solution  re¬ 
places  the  non-switch-optimal -solution  and  we  continue  with  the 
algorithm. 

(c)  Operation  2  for  &  =  1 :  A  switch  is  allowed  in  each  of  the  inter¬ 
vals  between  switches .  Since  there  is  no  further  improvement  in 
the  cost  function,  the  solution  for  l  *  0  is  the  optimal  solution. 
For  each  of  the  intervals  we  print  the  control  vector  u.  Note 
that  the  control  vector  u  is  printed  in  the  same  order  as  the 
capacities  and  the  columns  of  B. 

From  the  optimal  solution  and  from  matrix  B  we  conclude  that  in 
the  first  interval  Xj  ■  -2,  x^  ■  -2  and  x^  -  -3;  in  the  second  inter¬ 
val  Xj  -  0,  x2  *  -2  and  x^  ■  -3;  and  in  the  third  interval  x^  *  0, 
x2  *  -2  and  x^  *  0  (see  Figure  4.20). 


Figure  4.20:  State  trajectories  for  Example  4.1 


□  Example  4.1 


Figure  4.21:  Network  for  Example  4.2 
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The  results  of  the  program  are: 


THE  CAPACITY  VECTOR  IS: 

1 . 000 - i  .^OO  4.-0-QO  -  3^-003  - 2.  00a  5*0^00 — 


THE  8  MATRIX  IS: 

-1. 3-1  .0  0.0  0.0  0.0  0.0  0.0  0.0 

0.0  0*0- 1.0-1. 0  0.0  0.0  0.0  0.0 


0.0  0.0  0.0  0. 0-1.0  1.0  0.0  0.0 
0.0  0.0  0.0  0.0  0.0- 1.0- 1.0  0.0 
1.0  0.0  1.0  0.0  0.0  0.0  0. 0-1.0 

THE  INITIAL  CONDI T I UN  IS5 
1.00  4.00  21.00  15.00  20.00 

TAU<  1 ) =0.1 OOOOOOE  01 

U{  11=0.0 

U<  2 1=0 . 1 OOOOOOE  01 

■  J<  31=0.0 

HI  A ) =0. 300000 OE  01 

U<  5} =O.AOOOOOOK  01 

U(  61=0.0 

U  €  71 =O.SOOOOOOE  01 


U<  81 =0.A000000E  01 


tAUt  2  )=0  *333332  If-  00 

'J(  1 )=0.0 

U(  2  >-0.0 

!J{  3)  =0,0 

l(  A)  -0.3000000E  01 

J<  5)  =  0.  6000 00 0E  01 

*M  6) *0.0 

>><  7 )=o.*?ooooo it:  oi 

Ul  d  >  =0 • 4  00000 0£  01 

TAU(  3) =0.1 6666A7E  01 

>i(  1)=0.0 

'll  21=0,0 

H<  3)  =0  »  0 

il(  A  >  =  0.0 

U  (  5>=0. *0000005  01 

U(  6  1=0.0 

’J<  7)  =0.5000  00  OF  01 

IJC  8)=O.4000000E  01 


TACK  A)  =  0.5000000E  00 

uc  n=o.o 

U<  2  )  =0 . 0 

UC  31*0.0 

UC  4 >=0.0 

u(  5>=o.frooooooe  oi 

UC  6) =0.0 

UC  7)  =0.0 

UC  a> =0.40000026  Oi 

TAUC  6  }=0.if>00000E  01 

uc  n=o. o 

U<  2) =0.0 

UC  31=0.0 

UC  4) =0.0 

UC  51=0.0 

UC  6) =0. 0 

U C  71*0.0 

U(  8 1 *0 .4 00000 OE  Oi 
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Froa  the  optimal  solution  and  aatrix  B  we  conclude  that  (see  Figure  4.22): 

•  •  •  •  • 

first  interval:  Xj  •  -1,  x^  •  -3,  x3  *  -6,  x4  ■  -5*  xg  *  -4, 

second  interval:  Xj  •  0,  x^  *  -3,  x^  •  -6,  x^  *  -5,  Xj  *  -4, 

e  e  e  e  e 

third  interval:  x^  ■  0,  x^  a  0,  Xj  *  -6,  x4  *  -5,  xg  -4, 

•  •  e  •  * 

fourth  interval:  x^  *  0,  x^  •  0,  x^  «  -6,  x4  *  0,  xg  *  -4, 


□  Example  4.2 


Refer  to  the  network  appearing  in  Figure  4.23  with  the  initial  con¬ 
dition  x  =  (7,3,1,5,4,4,3)T. 

-o 


2 


-  8 

The  results  of  the  program  are: 

THE  CAPACITY  VECTOR  I  S  : 

3*000  A. 000  2*000  i .000  3.000  1*000  1.000  6.000  8.000 

a. ooo  2.000  e.ooc 

THE  e  MATRIX  is: 

-1.0— 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0 

0.0  0.0-1. 0-1.0  1.0  O.C  O.C  0.0  0.0  O.C  0.0  0.0 

0.0  0.0  0.0  O.O-l.O-l.C  1.0  C •  0  0.0  0.0  0.0  0.0 

O.C  0.0  0.0  0.0  0.0  l.C-l.C-1.0  0.0  0.0  0.0  0.0 

0.0  0.0  0.0  0.0  0.0  C.O  0.0  0. 0-1.0  1.0  0.0  0.0 

0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0-1. 0-1.0  0.0 

1  .0  0.0  1  .0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0— 1.0 

THE  INITIAL  CONDITION  IS! 


NC  SWITCHS  IN  THE  OPTIMAL  SOLUTION 
_ TA  UI  IJL?  Q,1  QjO  OO  OOe  41 - 

u(  » >*a.29999vae  oi 

UI  2IS0.A OOOOOAE  01 

U{  J)=3.  l99999dfE  01 

»J<  A»=O.IOOOOOOC  01 

U(  51=0.0 

U(  61 =0. 1 OQOOOOC  01 

U<  71*0.0 

<i(  si  =0. '‘■ooooooe  oi 

U(  91 =  0.80000961:  01 

UC  101  =  3. *00000*fc  01 

util)*- .3U146V7E- 05 

U  (  1  2  ) *0  »  79 V9999E  01 


In  this  work,  a  library  computation  program  for  solving  linear  pro¬ 
grams  was  used.  This  program  is  not  very  appropriate  for  solving  our  prob¬ 
lem.  For  small  networks  with  a  small  number  of  switches  in  the  optimal 
solution,  only  a  few  seconds  of  CPU  were  required  for  obtaining  the  optimal 
solution.  But  for  larger  networks  with  a  large  number  of  switches  in  the 
optimal  solution,  tens  and  even  hundreds  of  seconds  of  CPU  were  required. 
For  this  reason  larger  networks  could  not  be  completely  solved.  Further 
research  is  needed  for  obtaining  an  efficient  linear  program  that  will  per¬ 
mit  the  solution  of  large  networks. 


Open- Loop  Solutions  for  the  Dynamic  Routing  Problem 
by  Existing  Nonconvex  Quadratic  Programming  Methods 


5.1  Introduction 

In  this  chapter  we  shall  propose  other  approaches  for  solving  the 
open-loop  dynamic  routing  problem. 

In  Chapter  3,  the  open- loop  dynamic  routing  problem  was  presented 
as  a  nonconvex  quadratic  program  (Problem  3.2),  and  it  was  proved  that 
this  problem  is  not  concave,  quasi-concave,  quasi-convex  and  pseudo-convex 
(see  Section  3.3.3).  In  light  of  this  fact,  it  is  possible  to  solve  Problem 

3.2  by  using  existing  methods  for  nonconvex  quadratic  programming. 

In  this  chapter  we  will  indicate  the  most  appropriate  methods  with 
their  advantages  and  disadvantages. 

All  the  methods  mentioned  in  this  chapter  and  others  related  to 
nonconvex  quadratic  programming  are  over  viewed  in  Appendix  B. 

5.2  Ritter's  Algonthm 

This  algorithm  finds  a  global  minimum  of  a  nonconvex  quadratic  pro¬ 
gram  (see  Method  (19)  in  Appendix  B) .  It  is  a  finite  and  relatively  easily 
implementable  algorithm  for  finding  local  minima.  But,  in  spite  of  Ritter's 
convergence  proof  m  a  finite  number  of  steps  to  a  global  minimum,  Zwart 
found  an  example  that  can  not  be  solved  in  a  finite  number  of  steps.  There¬ 
fore,  in  order  to  solve  Problem  3.2  using  Ritter's  algorithm,  we  must  test 
convergence  for  this  specific  case. 
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5.3  The  Cabot-Francis  Algorithm 

Cabot  and  Francis  developed  an  algorithm  for  solving  the  nonconvex 
quadratic  programming  for  Thecases  inwhich  It  is  previously  known  that 
some  extreme  point  of  the  constraint  region  is  an  optimal  point.  Therefore, 
in  order  to  use  this  algorithm  for  the  solution  of  Problem  3.2,  it  must  be 
proved  that  the  above  mentioned  condition  holds.  For  this  purpose  we  need 
the  following  lemma: 

Lemma  5 . 1 

For  single  destination  networks  with  unity  weightings  in  the  cost 
functional,  there  exists  an  optimal  solution  to  the  dynamic  routing  problem, 
such  that: 

(a)  For  every  time  interval  [v  ^,tj),  i  e F,  the  optimal  control  is 
an  extreme  point  of 


Du* 

-l 

1?  * 

(5.1a) 

• 

u 

-l 

>  0  , 

(5.1b) 

•  T  * 
b  u. 

-j-i 

_0  , 

V  K  cl 

J  i 

(5.1c) 

a  o  , 

Vxf  E  B. 

1  1 

(5. Id) 

(5.1) 


(b)  The  number  of  switches  is  less  than  or  equal  to  n-1,  i.e., 

f  *:  n. 


Proof 

From  (2.20)  and  Proposition  2.2,  it  is  known  that  the  optimal  con 
trol  on  the  time  interval  ieF  is  obtained  from  the  linear 


program: 


— i  .  _  / _  ...  :  .  _ A 

AD-A0G5  109  MASSACHUSETTS  INST  OF  TECH  CAMBRIDGE  LAB  FOR  INFORMA— ETC  F/G  12/8 
OPEN-LOOP  SOLUTIONS  FOR  THE  DYNAMIC  ROUTING  PROBLEM* ( U I 
MAY  SO  S  SHATS*  A  SEGAL  N00014-7S*C-1183 

UNCLASSIFIED  LIDS-R-992  NL 
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subject  to 


T 

min  X  Bu^  , 

« 

Ui  eU  , 

i .  <  0  ,  V x.  e  I.  , 
x.  s  0  ,  Vx.  £  B.  . 

X»  X»  i  ^ 


C5.2a) 


(5.2b) 


From  the  linear  program  above  follows  that  there  exists  an  extreme  point  of 
(5.2b)  that  is  an  optimal  solution  to  (5.2).  Hence,  from  (2.5)  and  (2.7), 
follows  that  this  optimal  solution  is  an  extreme  point  of  (5.1). 


In  order  to  obtain  an  upper  bound  for  the  minimum  number  of  switches 
that' can  be  obtained  in  an  optimal  solution  satisfying  (5.1),  we  shall  look 
for  an  optimal  trajectory  with  maximum  number  of  switches.  In  order  to  do 
so,  we  assume  that  all  components  of  xq  are  greater  than  zero.  Starting 
from  tQ  we  apply  a  control  that  satisfies  part  (a)  of  the  lemma.  From 
Property  2.2  we  know  that  it  is  optimal  to  continue  with  this  control  until 
a  set  of  state  variables  become  zero,  or  until  x*(t)  reaches  the  origin. 
For  finding  the  upper  bound  to  the  minimum  number  of  switches,  we  assume 
that  only  one  state  variable  becomes  zero.  From  this  moment  we  apply  a  new 
control  until  another  state  variable  becomes  zero.  If  we  continue  in  this 
way  we  will  arrive  at  the  origin  with  maximum  n-  1  switches.  Then: 


or 


f  -  1  <_  n  -  1 
f  <  n 


□  Lemma  5 . 1 


This  lemma  allows  us  to  prove  Theorem  5.1  and  using  it  to  show  that 
Problem  3.2  has  the  property  required  by  the  Cabot-Francis  method. 

The  constraint  region  of  Problem  3.2  ((3.24)  -  (3.27))  is  presented 
in  (5.3)  in  a  different  way  because  we  want  to  associate  to  each  constraint 
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of  the  type  (5. Id)  a  slack  variable  equal  to  zero  and  to  each  constraint  of 
the  type  (5.1c)  a  possibly  nonzero  slack  variable  (see  (3.26)  and  (5.3c)). 
Moreover,  the  constraints  >_  0  are  not  written  in  (5.3)  because  from 
(3.24)  and  (3.25)  follows  that  they  are  redundant. 


Theorem  5.1 

In  single  destination  networks  with  unity  weightings  in  the  cost 
functional,  there  exists  an  optimal  solution  to  Problem  3.2  at  an  extreme 
point  of  the  constraint  region: 


Dz.  -  Ct.  <  0 
-x  -  l  — 


^  i  e  F 


(5.3a) 

(5.3b) 


P  P 

x  +  I  B5i_kl+  l  ki  =  0  '  P  “  2'3’--"f 

0  i-1  1  j=2  J 


(5.3c) 


f 

x  +  7  Bz .  =  0  , 
-o  . L.  -1 


(5.3d) 


(5.3) 


Proof 

All  the  extreme  points  of  a  linear  constraint  region  are  basic 
solutions  of  the  system  of  equations  defining  it.  And  each  basia  solution 
is  an  extreme  point.  Therefore,  a  point  is  an  extreme  point  if  and  only 
xf  the  number  of  nonzero  variables  is  less  than  or  equal  to  the  number  of 
constraints  defining  the  region.  The  theorem  will  therefore  be  proved  if 
we  show  that  there  exists  an  optimal  solution,  satisfying  Lemma  5.1,  that 
has  a  number  of  nonzero  variables  in  (5.3)  less  than  or  equal  to  the 
number  of  constraints  defining  (5.3). 


We  denote  by  V  the  constraint  region  of  all  control  variables 


u1(u2,...,u£  (formed  by  f  regions  of  the  type  (5.1))  and  define: 

o  =  (the  number  of  equations  of  the  type  (5.1a)  and  (5.1b)  in  V)  , 

8  =  (the  number  of  nonzero  variables  (including  slack  variables) 
in  the  equations  of  the  type  (5.1a)  and  (5.1b)  in  V)  . 

By  using  these  definitions  we  shall  count  the  number  of  nonzero 
variables  and  the  numbers  of  equations  in  V  for  an  optimal  solution. 

The  number  of  equations  of  the  type  (5.1c)  in  V  is  nf. 

If  we  denote  by  r  the  number  of  equations  of  the  type  b‘u.  *  0, 

““jL  “1 

i e  F,  the  number  of  nonzero  variables  (not  including  the  variables  u?  , 
i e F)  in  the  equations  of  the  type  (5.1c)  are  nf  -  r.  Note  that  all 
those  variables  are  slack  variables. 


From  the  above  mentioned  calculations  it  is  clear  that  for  an  optimal 
solution  satisfying  Lemma  5.1,  the  number  of  equations  in  V  is  nf  +  a  and 
the  number  of  nonzero  variables  in  V  is  8  +  nf-r. 


Since  the  optimal  solution  in  question  is  an  extreme  point  of  V 
(Lemma  5.1),  the  number  of  nonzero  variables  is  less  than  or  equal  to  the 
number  of  equations.  Therefore 


number  of  nonzero 
variables  in  V 


} 


8  +  nf-r  <  nf  +  a 


number  of 
equations  in  V 


which  implies 

6  -  r  <_  a  .  (5.4) 

We  now  prove  the  theorem,  by  counting  the  number  of  nonzero  variables 
and  equations  in  (5.3)  for  an  optimal  solution  satisfying  Lemma  5.1. 
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The  number  of  equations  of  the  type  (5.3a)  and  (S.3b)  is  equal  to 
the  number  of  equations  of  the  type  (5.1a)  and  (S.lb)  (i.e.  a),  since  they 
are  the  same  equations  multiplied  by  t^  ,  ieF.  Likewise  there  are 
n(f-l)  equations  in  (5. 3c)  and  n  equations  in  (5.3d).  Therefore,  the 
number  of  equations  in  (5.3)  is  nf_+_a. 


In  (5.3a)  and  (5.3b)  there  are  g  +  f  nonzero  variables;  the  f 
time  periods  ieF  (these  are  nonzero  since  «  0  means  that  there 

is  no  time  interval  [t^  an<*  the  5  nonzero  variables  of  equations 

(5.1a)  and  (S.lb)  (that  appeaT  in  (5.3)  multiplied  by  the  respective  t^). 

In  (5.3c)  there  are  n(f-l)  -  r  additional  nonzero  variables,  that 

is  all  the  slack  variables,  apart  from  the  r  slack  variables  associated 

T  * 

with  the  constraints  b^  =  0,  ieF. 

From  the  above  mentioned  facts  it  is  clear  that  the  number  of  non¬ 
zero  variables  in  (5.3)  is  g  +  f  * n(f-l)  -  r , 

From  Lemma  5.1  and  (5.4)  it  is  known  that 

g  -  r  a 
and 

f  _  n  , 

therefore 


(number  of  nonzero  \ 
(variables  in  (5.3)/ 


g  +  f  +  n(f-l) 


r 


nf  *  a 


(number  of  \ 

'(.equations  in  (5,3);' 


Consequently,  an  optimal  solution  satisfying  Lemma  5.1  is  an  extreme  point 


of  (5.3). 


□  Theorem  5 . L 


In  Theorem  5.1  it  was  proved  that  there  exists  an  optimal  solution 
to  Problem  3.2  that  is  an  extreme  point  of  (5.3).  Therefore  the  Cabot- 
Francis  method  can  be  used  for  solving  the  problem.  That  is,  the  problem 


can  be  solved  by  ranking  the  extreme  points  of  (S.3)  (see  Method  (22)  in 
Appendix  B) . 

The  most  serious  limitation  of  this  method  is  that  all  the  repre¬ 
sentations  of  some  of  the  extreme  points  of  (5.3)  are  needed  and  as  a  result 
of  that  the  efficiency  of  the  algorithm  is  not  good  for  degenerate  problems 
(Problem  3.2  can  be  very  degenerate). 

5.4  Other  Methods 

Besides  the  previously  mentioned  methods,  there  are  other  methods 
for  solving  nonconvex  quadratic  programs  that  may  be  used  for  the  solution 
of  Problem  3,2.  Some  of  those  methods  are:  the  Majthag,  Whinston  and 
Coffman  method  for  finding  a  local  minimum  (see  Method  (15)  in  Appendix  B) 
and  the  Candler-Townsley  method  for  finding  a  global  minimum  (see  Method 
(16)  in  Appendix  B) .  The  limitations  of  these  methods  are  that,  in  the 
first  only  a  local  minimum  is  reached,  and  in  the  second  there  is  no  proof 
of  convergence. 

Other  methods,  relatively  new,  that  assure  convergence  to  a  global 
minimum  are:  the  generalized  polars  method  (see  Method  (17)  in  Appendix  B) 
and  the  Tone  method  (see  Method  (18)  in  Appendix  B) . 


Conclusions 


6.1  Discussion 

In  this  work  we  made  use  of  the  fact  that  the  optimal  solution  is  of 
the  bang-bang  type,  in  order  to  formulate  the  open-loop  routing  problem  as 
dependent  on  new  variables,  u^  and  t^,  ieF,  The  problem  obtained  was 
shown  to  be  a  cubic  optimization  problem  having  a  nonconvex  constraint 
region.  Since  this  type  of  problems  is  quite  difficult,  the  fact  that  all 
nonlinear  terms  include  a  term  of  the  form  u.t^  was  used  for  transforming 
the  problem  into  a  quadratic  program.  This  new  problem  was  shown  to  be  non¬ 
convex,  with  nonconcave,  nonquasi -concave,  nonquasi -convex  and  nonpseudo- 
convex  cost  function,  and  therefore  it  does  not  possess  any  special  property 
that  will  allow  an  easier  solution  than  a  general  nonconvex  quadratic  program 

Consequently  two  approaches  were  proposed  for  the  solution  of  the 
problem.  The  first  one  based  on  existing  methods  from  nonconvex  quadratic 
programming,  the  second  consisting  of  developing  a  special  algorithm.  This 
special  algorithm  was  developed  for  single  destination  networks  with  unity 
weightings  in  the  cost  functional. 

The  algorithm  reduces  the  nonconvex  quadratic  program  to  a  series  of 
linear  programs.  The  number  of  switches  in  the  optimal  solution  depends  on 
the  initial  condition.  In  some  cases  the  number  of  switches  may  be  small 
and  one  obtains  relatively  small  linear  programs.  For  other  initial  condi¬ 
tions,  specially  in  big  networks,  the  optimal  solution  may  have  many  switches 
(up  to  n-1),  and  then  the  solution  of  large  linear  programs  may  be  required 
in  the  last  steps  of  the  algorithm. 


In  this  work  we  only  test  the  algorithm,  using  a  general -purpose, 
linear  program  from  a  library  of  computation  programs,  thus  excluding  the 
possibility  of  testing  large  problems. 

Further  research  is  required  in  order  to  develop  special  linear 
programs  that  are  adequate  for  solving  large  problems. 

6.2  Suggestions  for  Further  Work 

Some  suggestions  for  further  research  are: 

(a)  The  use  of  a  special-purpose  linear  program,  instead  of  the  general 
one  used  here,  for  obtaining  a  more  efficient  algorithm  that  will 
permit  solving  of  large  networks. 

(b)  All  the  results  of  this  work  refer  to  networks  without  inputs,  but 
they  can  be  extended  to  networks  with  constant  inputs 

(c)  The  algorithm  proposed  in  Chapter  4  was  developed  for  single¬ 
destination  networks  with  unity  weightings  in  the  cost  functional; 
further  research  is  needed  for  extending  the  results  to  multi - 
destination  networks  with  priorities  in  the  cost  functional. 

(d)  An  interesting  direction  for  further  investigation  may  be  to  develop 
an  algorithm  based  on  Ritter's  method,  that  solves  the  open-loop 
dynamic  routing  problem.  This  approach  may  be  interesting  because 
Ritter's  method  can  be  applied  to  general  networks  with  priorities 
in  the  cost  functional. 
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Appendix  A 

Properties  of  Functions 

A.l  Convexity  and  Concavity  of  Functions 
Definition  A.l 

A  set  S  is  said  to  be  convex  if  and  only  if 

Vx.yeS  then  (l-A)x  +  AyeS  for  Ae  [0,1]  . 

□  Definition  A,1 

Definition  A. 2 

A  real -valued  function  <t>  is  said  to  be  convex  on  S  if: 

(a)  S  is  convex  , 

(b)  Vx.yeS  then  (1-A)<J>  (x)  +  \<t>  Cy)  >_  [  Cl-A)x  +  Ay]  for  Ae  [0,1] 

□  Definition  A, 2 


Definition  A. 3 

A  real -valued  function  <f>  is  said  to  be  concave  on  A  if 

(a)  S  is  convex  t 

(b)  Vx.yeS  then  (1-A)<|>(x)  +  A<Ky)  <_  $[  (l-A)x  +  Ay]  for  Ac  [0,1] 

□  Definition  A. 3 

These  properties  are  important  in  mathematical  programming  since 
for  a  convex  function  the  Kuhn-Tucker  necessary  conditions  for  minimum  are 
also  sufficient  and  every  local  minimum  is  a  global  minimum.  Likewise,  if 
the  function  $  is  concave  on  S  the  minimum  is  achieved  at  an  extreme 
point  of  S. 
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We  now  define  some  properties,  weaker  than  convexity  and  concavity, 
that  are  sufficient  for  the  above  statements  to  hold. 

A. 2  Quasi-Convexity  and  Quasi-Concavity  of  Functions 
Definition  A. 4 

A  real -valued  function  $  defined  on  a  convex  set  S  is  said  to  be 
quasi -convex  on  S  if  and  only  if  the  set 

{xeS/4>(x)  ±  w}  is  convex  for  all  w  e  R  . 

□  Definition  A. 4 


This  property  is  partially  responsible  for  the  theoretical  useful¬ 
ness  of  quasi-convex  functions  as  constraints  in  a  nonlinear  program. 
Although  these  functions  need  not  be  convex,  they  still  determine  a  convex 
constraint  set. 


An  example  of  a  quasi-convex  function  that  is  not  convex  on  the 
interval  [0,1]  is 


<l>00 


Cj  if  0ix<i 


c2  if  2  <  x  i  1 


where  c^  and  are  scalars,  and  c^  ^  c^. 


Definition  A. 5 

A  real-valued  function  <J>  defined  on  a  convex  set  S  is  said  to  be 
quasi-convex  on  S  if  and  only  if 

Vx,yeS,  <fr[(l-A)x  +  Ay]  max[<Kx)  ,<Ky)]  for  le  [0,1]  . 

□  Definition  A. 5 


For  differentiable  functions,  the  following  definition  is  also 


equivalent: 
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Definition  A. 6 

A  differentiable  function  $  is  quasi -convex  on  the  convex  set  S 
if  and  only  if 

Vx.yeS,  $(y)  <_  $(x)  implies  V$(x) (y-x)  <_  0  , 

O  Definition  A. 6 

The  concept  of  quasi-convex  functions  was  introduced  by  De  Finetti 
[A8],  and  subsequently  developed  by  Nikaido  [A13],  and  by  Arrow  and  Ent- 
hoven  [A1 ] . 

Definition  A. 7 

A  real -valued  function  <f>  defined  on  a  convex  set  S  is  said  to  be 
quasi -concave  on  S  if  and  only  if  the  set 

{xeS/$(x)  '±.  w}  is  convex  for  all  weR  . 

□  Definition  A. 7 


Definition  A. 8 

A  real-valued  function  <{>  defined  on  a  convex  set  S  is  said  to  be 
quasi -concave  on  S  if  and  only  if 

Vx, yeS,  <J»[(1-X)x  +  Ay]  >  min[<Kx),<Ky)]  for  Xe  [0,1]  . 

□  Definition  A, 8 


Definition  A. 9 

A  differentiable  function  <|>  is  quasi-concave  on  the  convex  set  S 
if  and  only  if 

Vx,yeS,  <Ky)  <P(x)  implies  Vtfi(y)  (x-y)  >0  . 


□  Definition  A. 9 


An  important  property  of  quasi-concave  functions  is  that  the  mini¬ 
mum  of  the  function  4  is  achieved  at  an  extreme  point  of  the  convex  set  S. 
This  fact  can  be  easily  proved  by  noting  that  a  quasi-concave  function 
takes  on  only  higher  values  on  convex  combinations  than  it  takes  on  at  the 
defining  point  with  the  smallest  function  value. 

A. 3  Pseudo-Convexity  of  Functions 
Definition  A. 10 

A  real-valued  differentiable  function  $  on  a  convex  set  S  is 
pseudo-convex  if  and  only  if 

Vx.yeS,  7<)>(x)(y-x)  >;  0  implies  d> Cy)  >_  Cx)  . 

□  Definition  A. 10 

Every  local  minimum  of  a  pseudo-convex  function  is  global,  and  the 
Kuhn-Tucker  necessary  conditions  are  also  sufficient  for  a  local  (global) 
minimum  in  a  nonlinear  program  whose  objective  function  is  pseudo-convex 
and  the  constraints  are  defined  by  quasiconvex  functions.  See  [A2] . 

This  concept  was  first  introduced  by  Tui  [A19]  and  later  developed 
by  Mangasarian  [A9] ,  [A10] .  See  also  Ponstein  [A14] . 

Pseudo -convexity  is  defined  by  infinitely  many  inequalities  which 
are  in  general  hard  to  verify.  For  particular  functions,  such  as  products, 
quadratic  or  cubic  functions,  more  practical  criteria  are  available.  (See 
[A7] ,  [A10] ,  [A15] ,  [A16] ,  [A17]  and  [A18]. 


It  can  easily  be  shown  that: 


A  differentiable  convex  function  is  pseudo-convex  and  a  pseudo- 
convex  function  on  a  convex  set  is  quasi-convex. 

□  Proposition  A.l 

For  further  reading  on  the  convexity  of  quadratic  functions  see 
[A4],  [A5],  [A6] ,  and  on  the  quasi-convexity  and  pseudo-convexity  of 
quadratic  functions  see  [A3],  [A7],  [All],  [A12]  and  [A18] . 
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Nonconvex  Quadratic  Programming 


B.l  Introduction 

B.1.1  Definition  of  a  Quadratic  Program 

A  quadratic  program  can  be  defined  as  follows: 

T  T.. 

mm  z  =  c  x  +  x  Dx 

sub j  ect  to 


and 


Ax  *  b 


x  >_  0 


(B.l) 


where  x  and  c  are  column  vectors  of  dimension  n,  A  is  an  m*n  matrix, 
is  a  column  vector  of  dimension  m,  and  D  is  a  symmetric  n*n  matrix. 


Note  that  all  the  different  kinds  of  quadratic  programming  problems 
can  be  written  in  the  form  (B.l)  by  using  the  following  transformations: 

1.  It  is  possible  to  find  the  maximum  of  z  instead  of  the  minimum  by 
using 

max  z  =  min  -z  . 

_ 2.  It  is  possible  to  transform  inequality  constraints  into  equality 

constraints  by  adding  slack  variables. 

3.  It  is  possible  to  transform  a  set  of  i  free  variables  x j ,  j  e J 
into  l  *  1  nonnegative  variables  by  using  the  transformations: 

X !  a  X .  -  X  ,  jeJ, 

1  JO*' 

where  x .  >  0 ,  j  e  J ,  and  x  >  0 . 
j  —  o  — 
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B.1.2  Synopsis  of  the  Appendix 

In  Section  B.2  we  present  a  review  of  the  methods  for  solving  the  con¬ 
vex  quadratic  program,  i.e.  the  methods  for  solving  problem  (B.l)  when  D  is 
positive  definite  or  positive  semidefinite  (see  [A4],  [AS]  and  [A6]). 

In  Section  B.3  we  present  a  review  of  some  of  the  existing  methods  for 
solving  nonconvex  quadratic  programs.  This  will  be  done  by  classifying  the 
methods  into  four  groups  according  to  the  properties  of  the  cost  function. 

The  four  groups  are: 

1.  Methods  for  quasi-convex  or  pseudo-convex  cost  function  (see  B.3.1). 

2.  Methods  for  quasi-concave  or  concave  cost  function  (see  B.3. 2). 

3.  General  methods  (z  without  special  properties)  (see  B.3. 3). 

4.  Special  methods  (see  B.3. 4). 

For  each  one  of  these  groups  we  will  review  some  of  the  existing  methods. 

Since  the  problem  to  be  solved  (Problem  3.2)  has  properties  that 
permit  us  to  solve  it  by  using  Ritter’s  algorithm  (see  Section  B.3. 3),  or  the 
Cabot-Francis  algorithm  (see  Section  B.3.4),  we  shall  explain  these  methods 
in  depth. 

B.2  Convex  Quadratic  Programming 

Kuhn-Tucker ' s  Theorem  implies  that  a  necessary  condition  for  solving 
(B.l)  is  (see  [B13],  p.  2~13)~: 

v  -  A^X  -  2Dx  =  c 

x.v.  =  0,  j  -  1,2, . . . ,n  •  . 

Ax  *  b  , 
x  >  0,  v  >  .0 


(B.2) 
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In  convex  quadratic  programs,  each  stationary  poi,nt  is  a  global 
minimum,  therefore  the  conditions  (B.2)  are  also  sufficient  optimality 
conditions.  Almost  all  methods  for  solving  convex  quadratic  programs  are 
based  on  the  solution  of  this  system  of  equations . 

Some  of  the  most  important  methods  for  solving  convex  quadratic 
programs  are: 

(1)  Wolfe’s  Algorithm  (see  [B38];  and  [B30],  pp.  224-242) 

This  algorithm  is  based  ’upon  Kuhn-Tucker  conditions  and  the  SIMPLEX 
method  for  solving  linear  programs.  A  solution  to  (B.2)  is  achieved  by 
'•SIMPLEX- like"  steps  with  the  additional  requirement  that  x^v^  *  0  in 
each  step.  This  condition  is  known  as  the  complementarity  condition. 

(2)  Beale's  Method  (see  [B2] ;  and  [B30] ,  pp.  242-258) 

Beale's  method  is  based  (like  Method  (1))  on  the  Kuhn-Tucker  condi¬ 
tions  and  the  SIMPLEX  method,  but  in  this  method  the  constraint  region 
changes  during  the  application  of  the  algorithm. 

(3)  Lemke's  Method  (see  [B21];  and  [B30],  pp.  258-273) 

This  method  is  based  on  the  Kuhn-Tucker  conditions.  Since  D  *  is 
needed,  it  can  be  used  only  when  D  is  positive  definite. 

(4)  Theil  -  van  de  Panne's  Method  (see  [B34]) 

This  method  cons i s~t s  ~oT~TTnd ing  an -trpt-jmal-soTution^t o  the  problem 
without  constraints  and  then  adding  the  constraints  step  by  step,  where  in 
each  step  the  optimum  cost  is  calculated.  D  *  is  needed  for  the  calcula¬ 
tions,  therefore  this  method  can  be  used  only  when  D  is  positive  definite 


(5)  Capacity  Methods  (see  [B13],  pp.  236-238;  and  [B16J) 

These  methods  are  for  problems  with  constraints  of  the  type  Dx  <_  C 
One  of  them  is  the  Houthakker  method  that  is  based  on  the  introduction  of  a 
parametric  constraint  to  the  original  constraints. 

As  for  Methods  (3)  and  (4),  D  must  be  positive  definite. 

(6)  "Principle  Pivot11  Method  (see  [B8]) 

The  principle  pivot  method  is  a  method  for  solving  the  linear  com¬ 
plementarity  problem 

w  =  q  +  Mz 

wTz  =  0  (B.3) 

w  >  0,  z  >_  0 

The  Kuhn-Tucker  conditions  are  a  special  case  of  (B.3);  therefore, 
it  is  possible  to  solve  (B.2)  by  using  methods  for  solving  (B.3). 

In  [B14] ,  Hefley  compared  Methods  (1)  and  (6),  and  his  conclusion 
was  that  the  principle  pivot  method  as  appears  in  [B8]  is  superior. 

For  other  methods  see  [BIO],  [B12],  [B15] ,  [B18] ,  [B26]  and  [B29] . 


B.3  Nonconvex  Quadratic  Programming 

B.3,1  Special  Methods  for  Solving  Quadratic  Programs 

with  Pseudo-Convex _or _Quasi -Convex Tost  Functions -  - 

From  Section  A. 3  it  is  known  that  the  Kuhn-Tucker  conditions  are 
necessary  and  sufficient  for  minimum,  when  the  cost  function  is  pseudo- 
convex.  Therefore,  it  is  possible  to  solve  the  nonconvex  quadratic 
programs  with  pseudo-convex  cost  functions  by  using  some  of  the  methods 
of  Section  B.2. 


We  now  present  a  theorem  relating  pseudo-convex  with  quasi-convex 
functions : 

Theorem  B, 1  (see  [A3]) 

If  the  quadratic  function  z  (see  (B.l))  is  not  convex  but  is  quasi 
convex  in  the  nonnegative  orthant;  then,  if  c  ^  0,  the  function  z  is 
pseudo-convex  in  the  nonnegative  orthant. 

□  Theorem  B . 1 

Likewise,  in  [A12]  it  was  proved  that  when  c  =  0,  z  is  pseudo-convex  in 
the  nonnegative  orthant  only  when  z  is  convex. 

From  Proposition  A.l,  from  Theorem  B.l  and  from  the  above  follows 
that  pseudo-convex  and  quasi-convex  functions  are  closely  related  and 
almost  all  the  quadratic  quasi-convex  functions  are  also  pseudo-convex. 
Therefore  the  same  methods  can  be  used  for  solving  both  types  of  problems . 
For  the  case  where  the  cost  function  is  quasi-convex  and  c  *  0,  a  method 
for  solving  pseudo-convex  problems  can  be  used  if  a  small  perturbation  on 
c  is  performed. 

Martos  in  [All],  and  Mylander  in  [B25]  reviewed  the  methods  for 
solving  the  pseudo-convex  quadratic  programming  problem  and  concluded  that 
the  best  methods  are: 

-  3eaTe-s-method — (see  Method  (2)) 

Ritter's  method  (see  Method  (19)  in  Section  B.3.3) 

-  (7)  Keller's  Method  (see  [B18]) 

This  is  a  variation  of  Method  (6). 

(8)  Mylander 's  Method  (see  [B26] ) 

This  method  is  strongly  related  to  Lemke's  method  (see  Method  (3)). 


B.3.2  Special .Methods _f or  Solving  Quadratic  Programs 
with  Concave  or  <$uasi -Concave  Cost  Functions 

In  Sections  A.l  and  A. 2  it  was  shown  that  if  the  cost  function  is 
concave  or  quasi -concave,  the  minimum  of  the  function  is  located  at  an 
extreme  point  of  the  linear  constraint  region.  All  methods  for  solving 
these  problems  are  based  on  this  property. 

Some  of  the  most  important  methods  for  solving  optimization  problems 
with  concave  cost  function  (not  necessarily  quadratic)  are: 

(9)  Tui's  Method  (see  [B37]) 

Tui  was  one  of  the  first  investigators  in  dealing  with  this  problem. 
In  his  work  he  proposes  a  method  based  on  the  extension  of  the  cost  function 
(which  is  defined  on  a  specific  region)  over  Rn.  In  each  step  a  new  con¬ 
straint  is  created  (Tui  cut)  or  the  cost  function  is  changed. 

(10)  Hu's  Method  (see  [B17]) 

This  method  is  based  on  basical  elements  from  Tui's  work  and  on  the 
SIMPLEX  method.  During  the  steps  of  the  algorithm  the  constraint  region 
changes  since  the  method  requires  the  addition  of  constraints. 

(11)  Zwart 1 s  Method  (see  [B39]) 

This  algorithm  is  similar  to  Tui's  algorithm.  A  proof  of  convergence 
is  provided  for  the  case  when  a  little  deviation  from  the  optimum  is  permit¬ 
ted  and,  in  spite  of  the  fact  that  there  is  no  convergence  proof  for  the 
case  when  a  little  deviation  is  not  permitted,  the  programs  run  in  [B39] 
converged  to  the  optimum. 

(12)  Cabot's  Method  (see  [B4]) 

The  method  is  based  upon  the  Tui  cuts  (see  [B37])  and  on  the  branch- 


and-bound  technique. 
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(13)  Konno 1 s  Method  (see  [B20]) 

This  method  is  only  for  quadratic  functions  and  is  based  on  Konno' s 
method  for  solving  bilinear  problems  (see  Method  (20)). 

A  method  for  solving  optimization  problems  with  quasi-concave  cost 
function  (not  necessarily  quadratic)  is: 

(14)  Majthay  and  Whinston's  Method  (see  [B22]) 

This  method  is  based  on  placing  additional  cuts  to  the  constraint 
region  and  uses  the  fact  that  the  optimal  solution  is  located  at  an  extreme 
point  of  the  original  constraint  region.  Therefore,  the  algorithm  uses  a 
marking  method  for  differentiating  the  original  constraints  from  those  added 
by  the  algorithm. 

It  is  also  possible  to  solve  quadratic  programs  with  quasi-concave  or 
concave  cost  functions  by  using  Cabot  and  Francis'  method  (see  Method  (22)), 
or  by  using  all  the  methods  appearing  in  Section  B.3.3. 

B.3,3.  General  Methods  for  Solving  Nonconvex  Quadratic  Programs 

If  the  quadratic  function  has  none  of  the  properties  mentioned  in 
Sections  B.3.1  and  B.3.2,  the  quadratic  program  can  be  solved  by  using  one 
of  the  general  methods  for  solving  nonconvex  quadratic  programs. 

- Si  nett  Problem  3.2  is  a  nonconvex  quadratic  program,  the  methods _ 

appearing  in  this  section  are  more  extensively  explained;  special  attention 
is  given  to  Ritter's  algorithm,  since  it.  is  the  basic  work  on  the  subject. 

(15)  Majthay,  Whinston  and  Coffman's  Method  (see  [B23]) 

This  method  is  almost  identical  to  Ritter's  method  for  finding  a  local 
minimum,  with  the  only  difference  that  Cottle  and  Mylander's  presentation  is 
used  instead  of  the  original  Ritter  presentation  (see  Method  (19)). 


Phase  3 


As  said  before,  in  this  phase  a  cut  is  added.  This  cut  excludes  from 

the  constraint  region  all  regions  in  which  the  cost  function  cannot  achieve 

a  smaller  value  than  the  value  achieved  in  Phase  2.  This  is  done  by  finding 

a  global  minimum  to  a  new  problem  constructed  from  data  of  the  final  step 

T 

of  Phase  2.  This  new  problem  has  only  a  capacity  constraint  e  x  <_  r  and 
can  be  easily  solved  for  fixed  t.  Then  the  solution  to  this  problem  is 
used  for  the  construction  of  the  cutting  plane. 

The  algorithm  was  developed  by  Ritter  (see  [B27]  and  [B28]),  but  in 

[B9]  Cottle  and  Mylander  present  it  in  a  more  elegant  and  understandable  way. 

In  [B28]  Ritter  presents  a  proof  of  convergence  to  the  global  opti¬ 
mum  in  a  finite  number  of  steps,  but  in  [B40]  Zwart  shows  a  counterexample 
and  explains  why  Ritter's  proof  is  not  always  tTue. 

B.3.4  Methods  for  Solving  Problems  with  Cost  Functions 
having  Special  Properties 

Some  of  the  most  important  methods  for  solving  special  quadratic 
programs  are. 

(20)  Methods  for  Bilinear  Programs  (see  [B19]) 

A  bilinear  program  is  a  special  case  of  quadratic  programming  and 
can  be  defined  as  follows :  “  ~ 

T  T  T 

min  z(x,y)  -  c^  +  ^y  +  x  Cy  , 

subject  to 

A^x  =  bj,  A 2y  *  b2  ,  (B.7) 

x  0  ,  Y  l_Q  • 

One  of  the  methods  for  solving  this  problem  is  Konno's  method  (see 
[B19]).  This  method  is  the  basis  for  Konno's  algorithm  for  solving  concave 
quadratic  programs  (see  Method  (13) . 
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(21)  Method  for  Solving  the  Quasi -Concave  Problem 


T  T 

min  2  ■  (c  x+a) (d  x+$) 

s.t.  Ax  =  b  , 

x  >_  0 


(B.8) 


T  T 

where  (c  x+a)  and  (d  x+3)  are  strictly  positive  for  all  feasible 
solutions . 

The  method  was  developed  by  Swarup  (see  [B31],  [B32]  and  [B33]) 


(22)  Cabot  and  Francis1  Method  (see  [B5]  and  [B24]) 

There  are  certain  nonconvex  quadratic  programs  that  have  an  optimal 
solution  in  one  of  the  extreme  points  of  the  convex  set  of  feasible  solu¬ 
tions  (like  those  with  a  concave  or  quasi-concave  performance  index).  These 
problems  can  be  solved  by  Cabot  and  Francis '  algorithm  for  ranking  the 
extreme  points. 

In  Chapter  5  it  was  proven  that  Problem  3.2  has  the  required 
property,  therefore  it  will  be  explained  in  detail. 

In  this  method  the  extreme  points  are  ranked  in  increasing  order 
of  a  linear  function  received  from  the  quadratic  problem  (this  is  done  by 
Murty's  algorithm).  In  each  step  an  upper  and  lower  bound  to  the  optimum 
are  calculated  for  determining  whether  the  end  condition  is  satisfied  or 
if  it  is  necessary  to  continue  with  the  algorithm!  -  - - 

Before  presenting  the  algorithm,  we  present  its  theoretical  back¬ 
ground  and  Murty's  algorithm. 

The  problem  to  be  solved  is: 


_ jtfl 


Problem  B,1  (denoted  by  PI) 
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T  T 

min  f  (x)  *  c  x  +  x  Dx 


where 


S  =  (x/Ax  =  b,  x  >_  0} 


knowing  that  it  has  characteristics  that  justify  the  extreme  point  assumption. 

□  Problem  B.l 


If  we  denote  by  d^  the  j-th  column  of  D,  and  let 


u.  *  min  d.x 
3  -3- 


(B.9) 


s.t«  x  e  S  ,  Vj  =  1, . . . ,n  j 
we  may  state  the  boundary  linear  programming  problem: 


Problem  B.2  (denoted  by  P2) 


I*  y  ip 

lin  g (x)  =  l  (c.+u.)x.  =  (c  +u  )x 
j=i  J  J  J 


□  Problem  B.2 


It  can  be  easily  proved  that: 


Proposition  B.2 


g(x)  <  f(x),  VxeS 


□  Proposition  B.2 


And  as  a  consequence  of  this  we  can  state  that: 


Corollary  B.l 


If  x  is  an  optimum  solution  to  P2,  then  f„  *  g(x )  is  a  lower 
— o  •  -o 

bound  on  f*  (f*  ■  f (x*)  is  the  optimal  value  of  f) . 


E  Corollary  B.l 


Ill  - 


We  can  also  easily  compute  the  upper  bound  =  f  (xq) .  With  this 
upper  and  lower  bound  we  conclude: 

Corollary  B.2 

Given  any  upper  bound  fu  on  f*,  denote  by  {x^}  the  set  of  all 
extreme  points  x^  of  P2,  such  that  gQ^)  <_  fu.  Then  PI  has  an  optimum 
solution  such  that  x*e{x^}. 

□  Corollary  B.2 

Murty's  Algorithm 

This  algorithm  was  developed  to  solve  the  fixed  charge  problem: 


where 

and 

with 


and  has  two  parts: 


min  C(x)  *  D(x)  + z(x) 
xeS 


S  =  {x/Ax  =  b,  x  :>  0}  , 


n  n 

z(x)  =  I  c.x.,  D(x)  =  l  d  (1-6  )  (B. 10) 

j=l  J  J  -  j=i  3  °*xj 


6  =  < 
o,x. 


1  if  x.  »  0 
3 


(0  if  x.  >  0  , 


1 .  An  algorithm  for  ranking  the  vertices  of  S  in  increasing  order 

- of — - - — - - 


2.  An  algorithm  for  the  fixed  charge  problem,  assuming  that  the 

vertices  of  S  can  be  ranked  in  increasing  order  of  the  variable 
costs  z(x). 

For  our  purpose  we  are  interested  in  ranking  the  vertices  of  S 


gOO 


J-l 


3  3  3 


> 


in  increasing  order  of 
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that  is  to  say,  in  ranking  the  vertices  of  S  in  increasing  order  of  a 
linear  function;  then,  only  the  first  part  of  Murty's  method  is  relevant 
to  us: 

An  Algorithm  for  Ranking  the  Vertices  of  S  in  Increasing  Order  of  g(x) 
(first  part  of  Murty's  algorithm} 

Consider  the  linear  program  P2. 

n 

min  g(x)  *  I  (c.+u.)  ,  s.t.  xeS  . 

-  j*l  3  3 

Assuming  that  a  finite  optimum  exists,  it  is  well  known  that  there  exists  a 
vertex  that  is  optimal.  Consider  any  basic  solution  B,  corresponding  to 
any  nonbasic  variable  i  B,  let: 

-B  A 

Cj  =  relative  cost  of  coefficient  of  the  nonbasic  variable  x.  correspond¬ 
ing  to  the  basis  B. 

B  A 

8°  =  the  value  with  which  x.  enters  the  basis. 

B  A 

T.  =  the  new  basic  feasible  solution. 

Then  from  the  SIMPLEX  algorithm  follows  that 

g(T®)  -  g(B)  +  e®5®  , 

B 

where  Tj  represents  the  adjacent  extreme  points  to  the  extreme  point 
represented  by  B. 

Let  F(B)  denote  the  set  of  all  adjacent  vertices  of  B  with  cost 
value  not  less  than  that  of  B;  i.e., 

F(B)  *  {(T®:  0),  Vj/x  t  B,  E®  »  0}  .  (B.12) 

J  w  w 

But,  if  B  is  a  degenerate  basic  feasible  Solution,  let  Vg  denote  the 
vertex  represented  by  it  and  let  B^,...,Br  be  all  the  basic  solutions 
representing  the  same  vertex  V  .  Then  we  should  replace  (B.12)  by 

D 


(B.13) 


F<V 


r  B  B 

0^TiP:°)>  Vj/x.  *B,  C.p>0} 

p=l  •*  J  P  j 


Proposition  B.3 

Every  feasible  basic  solution  can  be  reached  by  taking  a  path  with 
nondecreasing  cost  (with  g(x)  as  the  cost  function)  from  xQ  through  the 
vertices  of  S  (where  xq  is  the  optimal  solution  of  P2) . 

□  Proposition  B.3 


Denoting  by  any  specific  basic  solution,  and  taking 
{x  ,...,x.  as  the  path  with  nondecreasing  cost  mentioned  in  Proposition 
B.3,  we  can  state: 

Proposition  B.4 

Suppose  (xq,...,x^  ^  are  already  known.  Let  us  define 

Fk  *  £PV  "  {?o,,,,’-k-l}  ' 

Then  x^  is  the  minimal  cost  solution  in  F^. 

□  Proposition  B.4 


Based  on  the  preceding  facts,  we  can  state  the  Cabot-Francis  algor¬ 


ithm  as  follows  (given  in  fldWodiagram  .form}: 


For  a  better  understanding  of  the  algorithm  we  will  solve  now  a 


small  example: 

Example  B.l 

Find  the  minimum  of 

,  12.  12 

f (x)  =  Xx  -  x2  -  ^  +  XjX2  +  fx2  , 

subject  to 

2Xj  +  x2  <_  4 

-x1  +  3x2<3  (B.  14) 

xx  >  0.  x2  >_  0 

The  constraint  region  is  represented  in  the  following  figure: 


Figure  B.2:  Constraint  region  for  Example  B.l 


For  solving  the  problem  we  must  express  (B.14)  in  the  standard  form 
of  PI,  and  this  is  done  by  adding  two  slack  variables  (Vj  and  v2) .  The 
resulting  problem  is: 


-  lie:- 


Find  the  minimum  of 


f(x)  =  (l,-l,0,0)x  +  x1 


-1/2 

0 

0 

0 


1  0 
1/2  0 
0  0 
0  0 


s.t. 


where  x  =  (x^x  ,v  ,v2). 


2x,  +  x„ 

+  v. 

*  4 

1  2 

1 

x,  +  3x„ 

+  V, 

*  3 

1  2 

2 

x 

>  0 

(B.15) 


In  order  to  obtain  the  boundary  linear  programming  problem  we  must 
calculate  u^  for  j  =  1,...,4. 


As  said  before, 


u.  *  mm 


J 


in  dTx  ,  Vj  *  1, . . . ,4  . 


xeS 


Hence : 


u  =  min  -  :  x 
1  xeS  1 


u  =  mm  x  +  =  x- 
Z  xeS  1  *  l 


u  =  mm 
3  xeS 


u,  =  min  0 
4  xeS 


Uj  *  -1  » 


u2  =  0  , 


u3  =  °  , 


u4  =  0  . 


Therefore  P2  is: 


s.t. 


min  g(x)  *  -x. 


x  e  S 


From  Figure  B.2  follows  that  the  minimum  fox  g(x)  is  achieved  at: 


X 

-o 


Now  we  proceed  with  the  algorithm  proposed  in  Figure  B.l. 


-  117  - 


1.  We  calculate  f.  and  £  . 

i  u 


ft  ’  s(50> 


f  =  f(X  ) 

u  -o' 


Note  that  the  SIMPLEX  tableau  representing  the  present  situation  (the 
tableau  corresponding  to  the  optimal  solution  of  P2)  is: 


1 

0 

3/7 

-1/7 

0 

9/7 

0 

1 

1/7 

2/7 

0 

10/7 

0 

0 

1/7 

2/7 

-1 

10/7 

2,  From  the  tableau  follows  that  the  next-best  extreme  point  is  received 
by  pivoting  with  ,  The  resulting  tableau  is: 


1 

-1/3 

0 

0 

1/3 

0 

0 

1/3 

-l" 

therefore,  g(Xj)  =  -1, 


3,  -i  /  ^  U(5l)  /  fu)  . 


hence 


=  -1 


■(Xj)  =  -1/2 


4.  -1/2  < 


(f(xx)  <  fu) 


fu--!/2  . 


hence 


From  the  tableau  follows  that  the  next-best  extreme  point  is  received  by 
pivoting  v2  with  x2*  The  resulting  tableau  is: 


therefore 


Hence,  x 


,  g(x2)  =  0  . 

(g(x2)  >  y 

is  the  optimal  solution: 


t  * 

<VV 

=  (0,1) 

and 

f(x*) 

=  -1/2 

□  Example  B.l 
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(16)  Candler  and  Townsley.'s  Method  (see  .  [B6J) 

In  this  method  a  SIMPLEX-based  technique  is  used  for  finding  a  local 
minimum  to  (B.l).  At  this  point  the  cost  function  is  represented  as  a  func¬ 
tion  of  the  nonbasic  variables  and  the  partial  derivatives  of  the  cost  func¬ 
tion  with  respect  to  the  nonbasic  variables  is  calculated.  These  derivatives 
are  the  basis  for  the  calculation  of  a  cut  that  eliminates  from  consideration 
the  local  minimum,  and  other  points  with  greater  value  of  the  cost  function, 
but  does  not  eliminate  points  with  smaller  cost  function  values.  Then  the 
steps  of  the  algorithm  are  reapplied. 

There  exists  no  convergence  proof  for  this  method.  In  [Bll ]  the 
proposed  method  was  used  for  solving  a  simple  example. 

(17)  Generalized  Polars  Method  (see  [Bl]  and  [B3] 

This  relatively  new  approach  is  based  upon  the  construction  of  the 
generalized  polar  of  the  Kuhn-Tucker  polyhedron  associated  with  the  quadratic 
program  (B,l).  This  generalized  polar  is  a  convex  polyhedral  cone  whose 
interior  contains  no  feasible  solution  better  than  the  best  known  one.  This 
and  other  properties  are  used  to  construct  a  polytope  containing  the  feasible 
set  and  contained  in  the  polar  of  the  latter.  This  is  done  by  starting  with 
a  simplex  defined  by  a  subset  of  the  problem  constraints,  and  successively 
activating  other  problem  constraints,  until  the  resulting  polytope  is  con¬ 
tained  in  the  polar  of  the  feasible  set.  When  this  is  achieved,  the  best 
complementary  solution  found  in  the  process  is  optimal,  or  none  exists. 

(18)  Tone's  Method  (see  [B35]) 

This  method  is  based  on  the  solution  of  the  generalized  linear  com¬ 


plementarity  problem,  received  from  the  Kuhn-Tucker  conditions  generalized 
by  the  addition  of  a  new  variable. 
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The  solution  of  the  generalized  linear  complementarity  problem  is 
based  on  an  algorithm  for  finding  all  extremal  rays  of  a  polyhedral  convex 
cone  with  some  complementarity  conditions  (see  [B36]). 

The  convergence  to  a  global  minimum  is  guaranteed, 

(19)  Ritter's  Algorithm  (see  [B9] ,  [B27],  [B28]  and  [B40]) 

This  algorithm  is  one  of  the  most  important  works  on  the  subject. 
This  fact  led  us  to  study  it  in  depth. 

The  algorithm  is  iterative  and  each  iteration  is  composed  of  the 
following  three  phases  (the  phases  will  be  explained  later): 

Phase  1:  Determines  a  vector  that  satisfies  the  constraints  and 
expresses  the  objective  function  in  terms  of  the  independent  (nonbasic) 
variables . 

Phase  2:  Determines  a  local  minimum,  with  the  extreme  point  that 
has  been  determined  in  Phase  1  as  a  starting  point. 

Phase  3:  Constructs  a  cutting  plane  that  excludes  the  previously 
located  local  minimum,  without  excluding  the  global  minimum  if  the  latter 
has  not  been  found  yet , 

After  the  cutting  plane  has  been  placed.  Phase  1  is  reapplied  to 
the  augmented  problem.  Termination  can  occur  in  Phase  1,  if  no  feasible 
points  remain  after  placing  the  cutting  plane;  in  Phase  2,  with  an  indica¬ 
tion  that  the  objective  function  is  not  bounded  below;  or  in  Phase  3,  if 
a  weak  sufficiency  condition  for  a  global  minimum  is  satisfied,  or  with  an 
indication  that  is  not  bounded  below  on  the  constraint  set. 


The  problem  to  be  solved  is: 


subject  to 


T  IT 

min  f (x)  =  c  x  +  Dx  * 


Ax  >  b  1 


x  >  0 


S  , 


(B  *4) 


or,  by  writing  the  constraints  more  concisely: 

subject  to 


T  IT 

min  f (x)  =  c  x  ♦  Dx 


Gx  >_  h 

For  its  solution  we  use  the  following  results: 


S  . 


(B.5) 


Lemma  B.l 

For  the  set  of  all  h  such  that  the  Kuhn -Tucker  conditions  of 
(B.5)  have  a  solution,  f  can  be  regarded  as  a  function  of  h  and  in  this 
sense 

«T 

i  -v  > 

a  being  the  Lagrange  multipliers  that  are  greater  than  iero,  and  h,  G 
being  the  appropriate  parts  of  h  and  G  respectively. 

□  Lemma  B . J 

Theorem  B.2 

If  (x*,A*)  is  a  solution  to  the  Kuhn-Tucker  conditions  of  (B.5), 
then  x*  is  a  local  minimum  of  f  on  S  iff  f  is  convex  on  the  set 

L  s  {x  e  Rn  /Gx  *  h) 

□  Theorem  B  2 

Proposition  B.l 

T 

If  f(x*)  is  a  local  minimum  of  (B.5),  and  min{x  Dx/G  x  0}  ^0 
»  "  xeS  ”  °" 

then,  f(x*)  is  a  global  minimum  of  f  on  S. 


Where  Gq  is  the  part  of  G  satisfying  the  constraints  as  equalities,  i.e. 


□  Proposition  B.l 

Note  that  Proposition  B.l  is  the  previously  mentioned  weak  sufficiency  condi 
tion  for  a  global  minimum. 

Using  the  above  results  we  proceed  to  explain  the  three  phases  of 
the  algorithm: 

Phase  1 

We  transform  the  inequality  constraints  to  equality  constraints  by 
adding  slack  variables  (v) . 

In  order  to  find  a  feasible  solution  we  add  artificial  variables  r. 

l 

n 

and  by  using  the  SIMPLEX  method  we  minimize  £  r. .  Termination  can  occur 

i=l  1 

with  a  feasible  point  (all  r.  «  0) ,  or  with  an  indication  that  there  is  no 
feasible  point. 

After  finding  the  feasible  point,  we  write  the  problem  in  such  a  way 
that  f  depends  only  on  the  nonbasic  variables.  It  must  be  noted  that  afte 
doing  that,  it  may  happen  that  f  will  depend  not  only  on  elements  of  x, 
but  on  elements  of  v  as  well. 

Now  we  redefine  x  in  such  a  way  that  x  ■  0  is  a  feasible  point. 

Phase  2 

As  said  before,  in  this  phase  a  local  minimum  is  found. 

The  previous  application  of  Phase  1  guaranties  that  x  =  0  is  a 


feasible  point.  The  algorithm  works  by  maintaining  a  solution  to  the  Kuhn- 
Tucker  conditions,  which  yields  to  a  local  minimum  for  (B.5)  augmented  with 
the  capacity  constraint: 
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T 

ex  <  t 


e.  >  0 
1  — 


(B.6) 


e.  >  0  if  c.  <  0 

l  —  l  — 


With  the  parameter  t  =  0,  the  point  x  =  0  is  obviously  a  local 
minimum  for  the  augmented  problem,  since  no  other  point  is  feasible.  At 
this  point,  the  capacity  constraint  is  binding  and  has  a  positive  multiplier. 
Then,  the  capacity  constraint  is  relaxed  by  an  increase  of  the  parameter 
and  the  Kuhn-Tucker  conditions  for  the  augmented  problem  are  used  to  compute 
a  local  minimum  x  =  x(t)  for  this  prohlem. 

However,  there  are  conditions  under  which  it  is  impossible  to  in¬ 
crease  r  so  that  x(t)  remains  a  local  minimum  of  the  augmented  problem. 
Depending  on  how  the  situation  arises,  it  will  be  necessary  either  to  specify 
that  certain  nonnegative  constraints  (x^  >_  0)  be  treated  as  equality  con¬ 
straints  (x^  =0),  or  to  change  the  capacity  constraints. 

Lemma  B,1  shows  that  as  r  is  being  increased,  the  objective  func¬ 
tion  is  being  decreased  as  long  as  the  capacity  constraint  has  a  positive 
multiplier.  After  a  finite  number  of  iterations,  the  algorithm  either 
terminates  with  an  indication  that  the  objective  function  has  no  finite 
lower  bound  on  the  feasible  region,  or  the  multiplier  of  the  capacity  con¬ 
straint  reaches  the  value  zero.  When  this  occurs,  we  have  a  local  minimum 
to  the  problem  without  the  capacity  constraint,  but  we  may  still  have  the 
imposed  condition  that  some  of  the  variables  be  held  at  the  value  zero.  If 
there  are  no  such  restrictions,  the  point  at  hand  is  a  local  minimum  to 
(B.5)  and  Phase  3  is  to  be  executed.  However,  if  we  have  a  local  minimum 
to  the  problem  with  the  additional  condition  that  some  of  the  variables  be 
equal  to  zero,  then  a  new  capacity  constraint  is  to  replace  the  one  just 
dropped  and  Phase  2  must  be  reapplied. 


for  Solving  the  Open-Loop  Optimal  Dynamic  Routing  Problem 
for  Single  Destination  Networks  with  Unity  Weightings  in  the  Cost  Functional 

C.l  Introduction 

The  program  presented  in  this  appendix  is  an  implementation  of  the 
algorithm  proposed  in  Chapter  4.  It  permits  to  solve  the  open-loop  optimal 
dynamic  routing  problem  for  single-destination  networks  with  unity  weightings 
in  the  cost  functional,  by  using  linear  programs  only  (see  Section  4.3). 

The  proposed  computation  program  is  written  in  FORTRAN  and  makes  use 
of  the  program  LA01B,  from  the  HARWELL  Library,  for  solving  linear  programs. 

C.2  The  Program  Input 

M  -  number  of  links  in  the  network 
N  -  number  of  nodes  in  the  network 

B  -  bidimensional  array  defining  the  network  topology  (x  *  Bu) .  Its 
dimension  is  (N,M) . 

CA  -  M-dimension  vector  representing  the  capacity  of  the  links 
XP  -  initial  condition  vector  of  dimension  N 

All  this  data  is  read  according  to  a  FORMAT  appearing  in  the  first  part  of 
the  program  (see  Section  C.6). 

IA  -  maximal  number  of  rows  in  the  A-matrix,  i.e.  o  (see  Section  C.4). 


This  data  is  inserted  directly  into  the  program  (see  Section  C.6). 
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C.3  The  Program  Output 

The  program  prints  first  the  problem  parameters,  i.e.,  the  capacity 
vector,  the  B-matrix  and  the  initial  condition  vector  (see  Section  C.6). 

Then  it  prints  the  optimal  solution,  i.e.,  the  optimal  flow  in  each  link  (the 
vector  u) ,  and  the  time  period  this  control  is  applied.  It  is  also  possible 
to  obtain  error  messages  from  the  linear  programs  and  intermediate  results 
from  the  algorithm,  by  using  IPR  =  1  instead  of  0. 

C.4  Dimensioned  Variables 

A  -  Bidimens ional  array  defining  the  constraint  region  of  the  linear 

program  to  be  solved  (see  Section  C.5).  Its  dimension  is  (a, 3)  where 

a  =  2N(L' +1)  +  M(2L'+2)  +  2L'  +  1  , 

6  *  M(2L'+2)  +  L'N  +  3L'  +  3  , 

L'  being  the  maximal  value  L  may  achieve  (L  =  0,1, 3, 7, . . . ,L') , 
Remember  that  the  number  of  switches  in  Operation  2  is  2L+  1  (see 
Figure  4,16). 

BD  -  Vector  defining  the  constraint  region  of  the  linear  program  to  be 
solved.  Its  dimension  is  a. 

C  -  Vector  of  the  cost  weightings.  Its  dimension  is  3. 

IND  -  Vector  required  by  LA01B.  Its  dimension  is  y  where: 

y  =  ct  tN(L'+l)  +M(2L'+2)  . 

WK  -  Vector  required  by  LA01B.  Its  dimension  is  y(a+3). 

TAU  -  Vector  of  the  time  intervals  =  t^  -  t^_j,  i  *  l,...,2fc+2  where 
*21+2  *  tf*  *tS  dimension  is  2L*  +  2. 
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X  -  Vector  of  the  variables  that  are  to  be  optimized  in  the  linear  program. 
Its  dimension  is  0. 

Z  -  Vector  of  dimension  M(2L'+2). 

SX  -  Vector  of  dimension  L' . 

SP  -  Vector  of  dimension  L*. 

C,5  Explanation  of  the  Problem 

In  the  first  part  of  the  program  the  parameters  of  the  problem  to  be 
solved  are  read  and  written.  In  the  second  part  the  constraint  region 
AX  <,  BP  is  built,  where 

'  H>  .  . 

*"(!)  ' 

Note  that  the  constraints  over  the  double  line  are  constraints  of  the  type 

X 

a  x  b  (the  LA01B  program  adds  slack  variables  and  transforms  them  into 

T 

equality  constraints),  and  a  cost  function  C  X  is  built  where 

••(:)  ■ 

It  is  easy  to  verify  that  the  linear  program  defined  above  corresponds 
to  Operation  1  of  the  algorithm  (see  (4.22)). 

After  solving  this  linear  program,  the  non-switch-optimal-solution, 
the  final  time  and  the  cost  function  value  are  stored  in  memory. 

In  the  third  part  of  the  program,  thfe  constraint  region  and  cost  func¬ 
tion  corresponding  to  Operation  2  are  built  (see  (4.23)),  and  the  resulting 
linear  program  is  solved.  The  corresponding  A,  BD,  C  and  X  are: 
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BD 


|  N(L+1) 

|  M(2L+2) 
}  N(L+1) 


L+l 


X  = 


-21 


TAU(l) 


_  TAU (2L+2)  J 


M(2L+2) 


LN 


2L+2 


CT  = 


T 

0  , 


-1  BTAU(l)  ' 


T 

0 

»  u  » •  •  • » 


-1  BTAU(L+1)' 


J 


rtT*  SxO-SXCir  n  SX(1)'  -  SX  (2)  • 

_  »  2  »u>  y  • 


M(2L+2) 

>  • 

,0joT 


LN 


2L  2 


L+l 

n 


11 

where  SXD  =  £  xD. ,  SX(i)'  =  ^  x!(t.').  The  "prime"  indicates  that  the 

i=l  1  j-1  3  1 

respective  value  is  known  from  the  previous  step  of  the  algorithm.  As  said 

before,  the  constraints  over  the  double  line  are  inequality  constraints  of 
T  ^ 

the  type  a  x  b. 


In  the  last  part  of  the  program,  the  optimal  cost  received  from  the 
last  application  of  Operation  2  is  compared  with  the  stored  value.  If  there 
is  no  decrease  in  the  cost  value,  the  optimal  solution  (i.e.,  the  stored 
solution)  is  printed  and  the  program  stops;  but  if  there  is  a  decrease  in  the 
optimal  cost,  L  is  increased  (L  ■  2L+1),  and  the  last  solution  and  cost 
value  are  stored  instead  of  the  previously  stored  solutions,  and  Operation  2 
is  reapplied.  If  there  is  some  TAU(i)'  equal  to  zero,  the  program  elimin- 


onn 


29 


ates  the  appropriate  elements  from  A,  BD,  C  and  X,  and  in  this  way  further 
switches  in  this  time  interval  are  not  permitted. 

C.6  The  Program 
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